text:
export decimation, masterTrigger, currentFrame, ramWriterMode, connectADC, startADC, stopADC, samplesPerPeriod, periodsPerFrame,
     numSlowDACChan, setSlowDACLUT, enableSlowDAC, slowDACStepsPerRotation, samplesPerSlowDACStep, prepareSlowDAC,
     currentWP, slowDACInterpolation, numSlowADCChan, numLostStepsSlowADC, bufferSize, keepAliveReset, triggerMode,
     slowDACStepsPerFrame, enableDACLUT, ADCPerformanceData, RPPerformance, RPStatus, RPInfo, startPipelinedData, PerformanceData, numChan

struct ADCPerformanceData
  deltaRead::UInt64
  deltaSend::UInt64
end

struct RPStatus
  overwritten::Bool
  corrupted::Bool
  stepsLost::Bool
  adcEnabled::Bool
  dacEnabled::Bool
end

struct PerformanceData
  wpRead::UInt64
  adc::ADCPerformanceData
  dac::DACPerformanceData
  status::RPStatus
end

struct RPPerformance
  data::Vector{PerformanceData}
end

struct RPInfo
  performances::Vector{RPPerformance}
end

function RPInfo()
  return RPInfo([RPPerformance([])])
end

decimation(rp::RedPitaya) = query(rp,"RP:ADC:DECimation?", Int64)
function decimation(rp::RedPitaya, dec)
  rp.decimation = Int64(dec)
  send(rp, string("RP:ADC:DECimation ", rp.decimation))
end

numChan(rp::RedPitaya) = 2

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

function setSlowDACLUT(rp::RedPitaya, lut::Array)
  lutFloat32 = map(Float32, lut)
  send(rp, string("RP:ADC:SlowDACLUT"))
  @debug "Writing slow DAC LUT"
  write(rp.dataSocket, lutFloat32)
end

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

function enableSlowDAC(rp::RedPitaya, enable::Bool, numFrames::Int64=0,
            ffRampUpTime::Float64=0.4, ffRampUpFraction::Float64=0.8)
  enableI = Int32(enable)
  return query(rp, string("RP:ADC:SlowDACEnable ", enableI,
              ",", numFrames, ",", ffRampUpTime, ",", ffRampUpFraction), Int64)
end

function slowDACInterpolation(rp::RedPitaya, enable::Bool)
  enableI = Int32(enable)
  send(rp, string("RP:ADC:SlowDACInterpolation ", enableI))
end

numSlowADCChan(rp::RedPitaya) = query(rp,"RP:ADC:SlowADC?", Int64)
function numSlowADCChan(rp::RedPitaya, value)
  send(rp, string("RP:ADC:SlowADC ", Int64(value)))
end

numLostStepsSlowADC(rp::RedPitaya) = query(rp,"RP:ADC:SlowDACLostSteps?", Int64)

function samplesPerPeriod(rp::RedPitaya)
  return rp.samplesPerPeriod
end
function samplesPerPeriod(rp::RedPitaya, value)
  rp.samplesPerPeriod = value
  samplesPerSlowDACStep(rp, value)
end

function periodsPerFrame(rp::RedPitaya)
  return rp.periodsPerFrame
end
function periodsPerFrame(rp::RedPitaya, value)
  rp.periodsPerFrame = value
end

samplesPerSlowDACStep(rp::RedPitaya) = query(rp,"RP:ADC:SlowDAC:SamplesPerStep?", Int64)
function samplesPerSlowDACStep(rp::RedPitaya, value)
  send(rp, string("RP:ADC:SlowDAC:SamplesPerStep ", value))
end

slowDACStepsPerRotation(rp::RedPitaya) = query(rp,"RP:ADC:SlowDAC:StepsPerRotation?", Int64)
function slowDACStepsPerRotation(rp::RedPitaya, value)
  send(rp, string("RP:ADC:SlowDAC:StepsPerRotation ", value))
end

function prepareSlowDAC(rp::RedPitaya, samplesPerStep, stepsPerRotation, numOfChan)
  numSlowDACChan(rp, numOfChan)
  samplesPerSlowDACStep(rp, samplesPerStep)
  slowDACStepsPerRotation(rp, stepsPerRotation)
end

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

function currentFrame(rp::RedPitaya)
  return Int64(floor(currentWP(rp) / (rp.samplesPerPeriod * rp.periodsPerFrame)))
end

function currentPeriod(rp::RedPitaya)
  return Int64(floor(currentWP(rp) / (rp.samplesPerPeriod)))
end

currentWP(rp::RedPitaya) = query(rp,"RP:ADC:WP:CURRENT?", Int64)
bufferSize(rp::RedPitaya) = query(rp,"RP:ADC:BUFFER:SIZE?", Int64)

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

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

# "TRIGGERED" or "CONTINUOUS"
function ramWriterMode(rp::RedPitaya, mode::String)
  send(rp, string("RP:RamWriterMode ", mode))
end

# "INTERNAL" or "EXTERNAL"
function triggerMode(rp::RedPitaya, mode::String)
  send(rp, string("RP:Trigger:Mode ", mode))
end

connectADC(rp::RedPitaya) = send(rp, "RP:ADC:ACQCONNect")
function startADC(rp::RedPitaya)
  send(rp, "RP:ADC:ACQSTATUS ON")
end
stopADC(rp::RedPitaya) = send(rp, "RP:ADC:ACQSTATUS OFF")

wasOverwritten(rp::RedPitaya) = query(rp, "RP:STATus:OVERwritten?", Bool)
wasCorrupted(rp::RedPitaya) = query(rp, "RP:STATus:CORRupted?", Bool)

function serverStatus(rp::RedPitaya)
  send(rp, "RP:STATus?")
  return readServerStatus(rp)
end

function readServerStatus(rp::RedPitaya)
  statusRaw = read!(rp.dataSocket, Array{Int8}(undef, 1))[1]
  status = RPStatus((statusRaw & 1) != 0, # overwritten
   (statusRaw & (1 << 1)) != 0, # corrupted
   (statusRaw & (1 << 2)) != 0, # stepsLost
   (statusRaw & (1 << 3)) != 0, # adcEnabled
   (statusRaw & (1 << 4)) != 0) # dacEnabled
  return status
end

function performanceData(rp::RedPitaya)
  send(rp, "RP:PERF?")
  return readPerformanceData(rp)
end

function readPerformanceData(rp::RedPitaya)
  adc = readADCPerformanceData(rp)
  dac = readDACPerformanceData(rp)
  return (adc, dac)
end

function readADCPerformanceData(rp::RedPitaya)
  perf = read!(rp.dataSocket, Array{UInt64}(undef, 2))
  return ADCPerformanceData(perf[1], perf[2])
end

# Low level read. One has to take care that the numFrames are available
function readSamples_(rp::RedPitaya, reqWP, numSamples)
  command = string("RP:ADC:DATA? ",Int64(reqWP),",",Int64(numSamples))
  send(rp, command)

  @debug "read data ..."
  u = read!(rp.dataSocket, Array{Int16}(undef, 2 * Int64(numSamples)))
  @debug "read data!"
  return u
end

# Low level read, reads samples, error and perf. Values need to be already requested
function readSamplesChunk_(rp::RedPitaya, reqWP::Int64, numSamples::Int64, into=nothing)
  @debug "read samples chunk ..."
  if into === nothing
    into = Array{Int16}(undef, 2 * Int64(numSamples))
  end
  data = read!(rp.dataSocket, into)
  status = readServerStatus(rp)
  (adc, dac) = readPerformanceData(rp)
  @debug "read samples chunk ..."
  perf = PerformanceData(UInt64(reqWP), adc, dac, status)
  return (data, perf)
end

# Low level read, that includes performance and error data
function readDetailedSamples_(rp::RedPitaya, reqWP::Int64, numSamples::Int64)
  command = string("RP:ADC:DATA:DETAILED? ",Int64(reqWP),",",Int64(numSamples))
  send(rp, command)
  return readSamplesChunk_(rp, reqWP, numSamples)
end

function readSamplesIntermediate_(rp::RedPitaya, reqWP::Int64, numSamples::Int64)
  data = readSamples_(rp, reqWP, numSamples)
  status = serverStatus(rp)
  (adc, dac) = performanceData(rp)
  perf = PerformanceData(UInt64(reqWP), adc, dac, status)
  return (data, perf)
end

function readSamplesOld_(rp::RedPitaya, reqWP::Int64, numSamples::Int64)
  data = readSamples_(rp, reqWP, numSamples)
  overwritten = wasOverwritten(rp)
  corrupted = wasCorrupted(rp)
  status = RPStatus(overwritten, corrupted, false, true, true)
  (adc, dac) = performanceData(rp)
  perf = PerformanceData(UInt64(reqWP), adc, dac, status)
  return (data, perf)
end

function startPipelinedData(rp::RedPitaya, reqWP::Int64, numSamples::Int64, chunkSize::Int64)
  command = string("RP:ADC:DATA:PIPELINED? ", reqWP, ",", numSamples, ",", chunkSize)
  send(rp, command)
end

# Low level read. One has to take care that the numFrames are available
function readDataSlow_(rp::RedPitaya, startFrame, numFrames)
  numPeriods = rp.periodsPerFrame
  numChan = numSlowADCChan(rp)

  command = string("RP:ADC:SLOW:FRAMES:DATA ",Int64(startFrame),",",Int64(numFrames))
  #send(rp, command)

  @debug "read data ..."
  #u = read!(rp.dataSocket, Array{Float32}(undef, numChan * numFrames * numPeriods))
  @debug "read data!"
  #return reshape(u, numChan, numPeriods, numFrames)
  return zeros(Float32, numChan, numPeriods, numFrames)
end

max_stars_repo_name: Fypsilonn/RedPitayaDAQServer

max_issues_repo_name: Fypsilonn/RedPitayaDAQServer

max_forks_repo_name: Fypsilonn/RedPitayaDAQServer


text:
using PackageCompiler

import Pkg
Pkg.activate(".")

PackageCompiler.create_sysimage([:Deductive]; sysimage_path="sys_deductive.so", precompile_execution_file="precompile_execution.jl")

ext: jl


text:
module UnitRots
    include("Krav.jl")
    include("Emodes.jl")
    include("AuxRep.jl")
end

max_stars_repo_name: rurz/UnitRots

max_issues_repo_name: rurz/UnitRots

max_forks_repo_name: rurz/UnitRots


text:
using ChaosTools
using Test
using Distances: Cityblock, Euclidean

test_value = (val, vmin, vmax) -> @test vmin <= val <= vmax

println("\nTesting nonlinear timeseries analysis...")
@testset "Henon reconstruct" begin
    ds = Systems.henon()
    data = trajectory(ds, 100000)
    x = data[:, 1] # some "recorded" timeseries

    @testset "Dimension" begin
        τ = 1; D = 2
        R = reconstruct(x, D-1, τ)
        D2 = information_dim(R)
        test_value(D2, 1.1, 1.3)
    end
    @testset "Numerical Lyapunov" begin
        ks = 1:20
        @testset "meth = $meth" for meth in
            [FixedMassNeighborhood(1), FixedMassNeighborhood(4),
            FixedSizeNeighborhood(0.01)]
            @testset "distance = $di" for di in [Euclidean(), Cityblock()]
                for D in [2, 4]
                    R = reconstruct(x, D-1, 1)
                    E = numericallyapunov(R, ks,
                    refstates = 1:1000, distance=di, ntype=meth)
                    λ = linear_region(ks, E)[2]
                    test_value(λ, 0.3, 0.5)
                end
            end
        end
    end
end
@testset "Multidim Multitime R" begin
    ds = Systems.towel()
    data = trajectory(ds, 10000)

    taus = [2 3; 4 6; 6 8]
    data2 = data[:, 1:2]
    R = reconstruct(data2, 3, taus)

    ks = 1:20
    E = numericallyapunov(R, ks,
    refstates = 1:1000, distance=Cityblock(), ntype=FixedMassNeighborhood(1))
    λ = linear_region(ks, E)[2]
    test_value(λ, 0.3, 0.6)
end

@testset "Broomhead-King" begin
    using Random
    Random.seed!(1234);
    ds = Systems.gissinger(ones(3)) # standard initial condition

    data = trajectory(ds, 1000.0, dt = 0.05)
    x = data[1:end-1, 1] # "exactly" 20000 points
    s = x .+ 0.01rand(length(x))

    Ux, Σx = broomhead_king(x, 40)
    Us, Σs = broomhead_king(s, 40)
    R = reconstruct(x, 1, 30)

    for j in 10:40
        @test Σx[j] < Σs[j]
    end

    newcoords = Dataset(Us[:, 1], Us[:, 2])
    newcoordsclean = Dataset(Ux[:, 1], Ux[:, 2])

    Dnew = information_dim(newcoords)
    DR = information_dim(R)
    DC = information_dim(newcoordsclean)
    @test abs(Dnew - 0.1 - DR) < 0.2 # subtract 0.1 for "added dimensionality"
    @test abs(DC - DR) < 0.2
end


text:
export
    rotations

rotations(v::AbstractArray{T,1}) where {T} = VectorRotationsIterator(v)

struct VectorRotationsIterator{T}
    v::T
end

Base.start(::VectorRotationsIterator) = 0
Base.next(it::VectorRotationsIterator, state) = circshift(it.v, state), state - 1
Base.done(it::VectorRotationsIterator, state) = abs(state) >= length(it.v)

Base.eltype(it::VectorRotationsIterator) = Base.eltype(typeof(it))
Base.eltype(::Type{VectorRotationsIterator{T}}) where {T} = T

Base.length(it::VectorRotationsIterator) = length(it.v)


text:
reference_site = "http://www.bic.mni.mcgill.ca/ServicesAtlases/ICBM152Lin"
download = "http://packages.bic.mni.mcgill.ca/mni-models/icbm152/mni_icbm152_lin_nifti.zip"


text:
# AbstractDiscretizer
abstract type AbstractDiscretizer{N,D} end
    # N indicates the decoded, or natural type
    # D indicates the encoded, or discrete type

encoded_type(::AbstractDiscretizer{N,D}) where {N,D} = D
decoded_type(::AbstractDiscretizer{N,D}) where {N,D} = N

# DiscretizationAlgorithm
abstract type DiscretizationAlgorithm end

function get_discretization_counts(disc::AbstractDiscretizer, data::AbstractArray)
    counts = zeros(Int, nlabels(disc))
    for v in data
        counts[encode(disc, v)] += 1
    end
    counts
end

function get_discretization_counts(discs::Vector{T}, data::AbstractMatrix{N}) where {N, D, T<:AbstractDiscretizer{N, D}}
    nobservations, ndimensions = size(data)
    @assert length(discs) == ndimensions

    nbins = Tuple(nlabels.(discs))
    counts = zeros(Int, nbins)
    for i in 1:nobservations  # for each data point
        binaddress = zeros(D, ndimensions)
        for j in 1:ndimensions  # for each dimension
            binaddress[j] = encode(discs[j], data[i, j])
        end
        counts[binaddress...] += 1
    end

    counts
end

function get_histogram_plot_arrays(binedges::Vector{R}, counts::AbstractVector{I}) where {R<:Real, I<:Integer}
    n = length(binedges)
    n == length(counts)+1 || error("binedges must have exactly one more entry than counts!")
    arr_x = [binedges[1]; binedges]
    arr_y = [zero(I); counts; zero(I)]
    (arr_x, arr_y)
end
function get_histogram_plot_arrays(binedges::Vector{R}, pdfs::AbstractVector{F}) where {R<:Real, F<:Real}
    n = length(binedges)
    n == length(pdfs)+1 || error("binedges must have exactly one more entry than pdfs!")
    arr_x = [binedges[1]; binedges]
    arr_y = [zero(F); pdfs; zero(F)]
    (arr_x, arr_y)
end


text:
##
using QuasiMonteCarlo: LowDiscrepancySample, sample
using Plots

lb = [0.0, -0.5]  # Lower bounds
ub = [4.0, 3.5]   # Upper bounds
n = 100           # Number of samples

halton = LowDiscrepancySample([3, 5])  # Halton with prime number sequence
xtrain = sample(n, lb, ub, halton)

scatter(xtrain[1,:], xtrain[2,:])
##

# f(x) = cos.(x[1,:]).*sin.(x[2,:])
# ytrain = f(xtrain)

# #scatter3d(xtrain[1,:], xtrain[2,:], ytrain, marker_z=ytrain)
# #savefig("test.png")

max_stars_repo_name: krystophny/profit

max_issues_repo_name: krystophny/profit

max_forks_repo_name: krystophny/profit


text:
  # Function to be minimized
  f(x::Vector{Float64}) = x[1]^2 + x[2]^2

  # Minimizer by random search
  function randomsearch(f,ntrial)
    fbest = 1.e30
    x = Vector{Float64}(undef,2)
    xbest = Vector{Float64}(undef,2)
    xbest[1] = -10. + 20. * rand()
    xbest[2] = -10. + 20. * rand()
    for i in 1:ntrial
      x[1] = xbest[1] + 1.e-3*(-1.e0 + 2.e0 * rand())
      x[2] = xbest[2] + 1.e-3*(-1.e0 + 2.e0 * rand())
      fx = f(x)
      if fx < fbest
        fbest = fx
        xbest[1] = x[1]
        xbest[2] = x[2]
        println(i," New best point: ", x," f(x) = ", fx)
      end
    end
    println(" Best point found: ",xbest," f = ", fbest)
  end


text:
######### Stan batch program example  ###########

using Stan, Mamba

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  eightschools ="
  data {
    int<lower=0> J; // number of schools
    real y[J]; // estimated treatment effects
    real<lower=0> sigma[J]; // s.e. of effect estimates
  }
  parameters {
    real mu;
    real<lower=0> tau;
    real eta[J];
  }
  transformed parameters {
    real theta[J];
    for (j in 1:J)
      theta[j] <- mu + tau * eta[j];
  }
  model {
    eta ~ normal(0, 1);
    y ~ normal(theta, sigma);
  }
  "

  schools8data = [
    Dict("J" => 8,
      "y" => [28,  8, -3,  7, -1,  1, 18, 12],
      "sigma" => [15, 10, 16, 11,  9, 11, 10, 18],
      "tau" => 25
    )
  ]

  global stanmodel, rc, sim1
  stanmodel = Stanmodel(name="schools8", model=eightschools);
  rc, sim1 = stan(stanmodel, schools8data, ProjDir, CmdStanDir=CMDSTAN_HOME)

  if rc == 0
    nodesubset = ["lp__", "accept_stat__", "mu", "tau", "theta.1", "theta.2", "theta.3", "theta.4", "theta.5", "theta.6", "theta.7", "theta.8"]

    ## Subset Sampler Output
    sim = sim1[1:size(sim1, 1), nodesubset, 1:size(sim1, 3)]

    describe(sim)
    println()

    ## Brooks, Gelman and Rubin Convergence Diagnostic
    try
      gelmandiag(sim1, mpsrf=true, transform=true) |> display
    catch e
      #println(e)
      gelmandiag(sim, mpsrf=false, transform=true) |> display
    end

    ## Geweke Convergence Diagnostic
    gewekediag(sim) |> display

    ## Highest Posterior Density Intervals
    hpd(sim) |> display

    ## Cross-Correlations
    cor(sim) |> display

    ## Lag-Autocorrelations
    autocor(sim) |> display

    ## Deviance Information Criterion
    #dic(sim) |> display

    ## Plotting

    p = plot(sim, [:trace, :mean, :density, :autocor], legend=true);
    draw(p, nrow=4, ncol=4, filename="$(stanmodel.name)-summaryplot", fmt=:svg)
    draw(p, nrow=4, ncol=4, filename="$(stanmodel.name)-summaryplot", fmt=:pdf)
  end # rc == 0
end # cd


text:
## Table

function rechunk_together(left, right, lkey, rkey,
                          lselect=excludecols(left, lkey), rselect=excludecols(right, rkey); chunks=nworkers())
    # we will assume that right has to be aligned to left
    l = reindex(left, lkey, lselect)
    r = reindex(right, rkey, rselect)

    if has_overlaps(l.domains, closed=true)
        l = rechunk(l, chunks=chunks)
    elseif !issorted(l.domains, by=first)
        perm = sortperm(l.domains, by=first)
        l = fromchunks(l.chunks[perm], domains=l.domains[perm])
    end

    splitters = map(last, l.domains)

    r = rechunk(r, rkey, rselect,
                splitters=splitters[1:end-1],
                chunks_presorted=true,
                #affinities=map(x->first(Dagger.affinity(x))[1].pid, l.chunks),
               )
    l, r
end

function Base.join(f, left::DDataset, right::DDataset;
                   how=:inner,
                   lkey=pkeynames(left), rkey=pkeynames(right),
                   lselect=left isa DNDSparse ? valuenames(left) : excludecols(left, lkey),
                   rselect=right isa DNDSparse ? valuenames(right) : excludecols(right, rkey),
                   broadcast=nothing,
                   chunks=nworkers(),
                   kwargs...)
    cl = compute(left)
    cr = compute(right)

    if broadcast === :left
        error("only broadcast = :right is supported at the moment")
    elseif broadcast === :right
        if !(how in [:inner, :left, :anti])
            error("Can only do inner, left, or anti join with broadcast = :right")
        end
        if length(cr.chunks) == 1
            right_ser = cr.chunks[1]
        else
            right_ser = collect(cr)
        end
        ps = map(cl.chunks) do c
            a = Dagger.affinity(c)
            first.(a[1:min(length(a), 1)])
        end
        tasks = [delayed(identity)(right_ser) for p in ps]
        for (t, p) in zip(tasks, ps)
            t.affinity = Nullable(Pair.(p, 1))
        end
        r = fromchunks(tasks)
        l = cl
    else
        l, r = rechunk_together(left, right, lkey, rkey, lselect, rselect,
                                chunks=chunks)
    end

    i = 1
    j = 1

    cs = []
    rempty = delayed(empty!∘copy)(r.chunks[end])
    lempty = delayed(empty!∘copy)(l.chunks[end])

    function joinchunks(x, y)
        res = if broadcast !== nothing
            join(f, x, y, lkey=lkey, rkey=rkey,
                 lselect=lselect, rselect=rselect,
                 how=how; kwargs...)
        else
            join(f, x, y, how=how; kwargs...)
        end
    end

    while i <= length(l.chunks) && j <= length(r.chunks)
        c1 = l.chunks[i]
        c2 = r.chunks[j]

        pkey_tuple(nt,pkey) = tuple( ( getproperty(nt,n) for n in pkey )... )
        pkey_domain(d,pkey) =
            Interval(pkey_tuple(first(d),pkey),pkey_tuple(last(d),pkey))
        d1 = pkey_domain(l.domains[i].interval, lkey)
        d2 = pkey_domain(r.domains[j].interval, rkey)
        if hasoverlap(d1, d2)
            i += 1
            j += 1
        elseif last(d1) < last(d2)
            # this means there's no corresponding chunk on the right
            c2 = rempty
            i += 1
        else
            # this means there's no corresponding chunk on the left
            c1 = lempty
            j += 1
        end
        c = delayed(joinchunks)(c1,c2)
        push!(cs, c)
    end

    while i <= length(l.chunks)
        push!(cs, delayed(joinchunks)(l.chunks[i], rempty))
        i+=1
    end

    while j <= length(r.chunks)
        push!(cs, delayed(joinchunks)(lempty, r.chunks[j]))
        j += 1
    end

    fromchunks(cs)
end

function Base.join(left::DDataset, right; how=:inner, kwargs...)
    f = how === :anti ? ((x,y)->x) : IndexedTables.concat_tup
    join(f, left, right; how=how, kwargs...)
end

function groupjoin(f, left::DDataset, right; how=:inner, kwargs...)
    join(f, left, right; how=how, group=true, kwargs...)
end

function groupjoin(left::DDataset, right; how=:inner, kwargs...)
    join(left, right; how=how, group=true, kwargs...)
end

function join(f, left::DDataset, right::IndexedTables.Dataset; how=:inner, kwargs...)
    if how in [:inner, :left, :anti]
        join(f, left, distribute(right, 1), broadcast=:right, how=how; kwargs...)
    else
        join(f, left, distribute(right, 1), how=how; kwargs...)
    end
end

## NDSparse join

"""
    naturaljoin(left::DNDSparse, right::DNDSparse, [op])

Returns a new `DNDSparse` containing only rows where the indices are present both in
`left` AND `right` tables. The data columns are concatenated.
"""
function naturaljoin(left::DNDSparse{I,D1}, right::DNDSparse{J,D2}) where {I,J,D1,D2}
    naturaljoin(IndexedTables.concat_tup, left, right)
end

"""
    naturaljoin(op, left::DNDSparse, right::DNDSparse, ascolumns=false)

Returns a new `DNDSparse` containing only rows where the indices are present both in
`left` AND `right` tables. The data columns are concatenated. The data of the matching
rows from `left` and `right` are combined using `op`. If `op` returns a tuple or
NamedTuple, and `ascolumns` is set to true, the output table will contain the tuple
elements as separate data columns instead as a single column of resultant tuples.
"""
function naturaljoin(op, left::DNDSparse, right::DNDSparse)
    out_domains = Any[]
    out_chunks = Any[]

    # if the output data type is a tuple and `columns` arg is true,
    # we want the output to be a Columns rather than an array of tuples
    for i in 1:length(left.chunks)
        lchunk = left.chunks[i]
        subdomain = left.domains[i]
        lbrect = subdomain.boundingrect
        # for each chunk in `left`
        # find all the overlapping chunks from `right`
        overlapping = map(boundingrect.(right.domains)) do rbrect
            boxhasoverlap(lbrect, rbrect)
        end
        overlapping_chunks = right.chunks[overlapping]

        # each overlapping chunk from `right` should be joined
        # with the chunk `lchunk`
        joined_chunks = map(overlapping_chunks) do r
            delayed(naturaljoin)(op, lchunk, r)
        end
        append!(out_chunks, joined_chunks)

        overlapping_domains = map(r->intersect(subdomain, r),
                                     right.domains[overlapping])

        append!(out_domains, overlapping_domains)
    end
    return fromchunks(out_chunks, domains = out_domains)
end

Base.map(f, x::DNDSparse{I}, y::DNDSparse{I}) where {I} = naturaljoin(x, y, f)

# left join

"""
    leftjoin(left::DNDSparse, right::DNDSparse, [op::Function])

Keeps only rows with indices in `left`. If rows of the same index are
present in `right`, then they are combined using `op`. `op` by default
picks the value from `right`.
"""
function leftjoin(op, left::DNDSparse{K,V}, right::DNDSparse,
                  joinwhen = boxhasoverlap,
                  chunkjoin = leftjoin) where {K,V}

    out_chunks = Any[]

    for i in 1:length(left.chunks)
        lchunk = left.chunks[i]
        subdomain = left.domains[i]
        lbrect = subdomain.boundingrect
        # for each chunk in `left`
        # find all the overlapping chunks from `right`
        overlapping = map(rbrect -> joinwhen(lbrect, rbrect),
                          boundingrect.(right.domains))
        overlapping_chunks = right.chunks[overlapping]
        if !isempty(overlapping_chunks)
            push!(out_chunks, delayed(chunkjoin)(op, lchunk, treereduce(delayed(_merge), overlapping_chunks)))
        else
            emptyop = delayed() do op, t
                empty = NDSparse(similar(keys(t), 0), similar(values(t), 0))
                chunkjoin(op, t, empty)
            end
            push!(out_chunks, emptyop(op, lchunk))
        end
    end

    cache_thunks(DNDSparse{K,V}(left.domains, out_chunks))
end

leftjoin(left::DNDSparse, right::DNDSparse) = leftjoin(IndexedTables.concat_tup, left, right)

function asofpred(lbrect, rbrect)
    allbutlast(x::Interval) = Interval(first(x)[1:end-1],
                                       last(x)[1:end-1])
    all(boxhasoverlap(lbrect, rbrect)) ||
    (all(boxhasoverlap(allbutlast(lbrect),
                       allbutlast(rbrect))) &&
     !isless(last(lbrect), first(rbrect)))
end

function asofjoin(left::DNDSparse, right::DNDSparse)
    leftjoin(IndexedTables.right, left, right, asofpred, (op, x,y)->asofjoin(x,y))
end

function merge(left::DNDSparse{I1,D1}, right::DNDSparse{I2,D2}; agg=IndexedTables.right) where {I1,I2,D1,D2}

    I = promote_type(I1, I2)        # output index type
    D = promote_type(D1, D2)        # output data type

    t = DNDSparse{I,D}(vcat(left.domains, right.domains),
                    vcat(left.chunks, right.chunks))

    overlap_merge(x, y) = merge(x, y, agg=agg)
    if has_overlaps(t.domains)
        t = rechunk(t,
                    merge=(x...)->_merge(overlap_merge, x...),
                    closed=agg!==nothing,
                    sortchunks=false)
    end

    return cache_thunks(t)
end

function merge(left::DIndexedTable, right::DIndexedTable; chunks=nworkers())
    l, r = rechunk_together(left, right, pkeynames(left), pkeynames(right); chunks=chunks)
    fromchunks(delayedmap(merge, l.chunks, r.chunks))
end

function subbox(i::Interval, idx)
    Interval(i.first[idx], i.last[idx])
end

function bcast_narrow_space(d, idxs, fst, lst)
    intv = Interval(
        tuplesetindex(d.interval.first, fst, idxs),
        tuplesetindex(d.interval.last, lst, idxs)
    )

    box = Interval(
        tuplesetindex(d.boundingrect.first, fst, idxs),
        tuplesetindex(d.boundingrect.last, lst, idxs)
    )
    IndexSpace(intv, box, Nullable{Int}())
end

function broadcast(f, A::DNDSparse, B::DNDSparse; dimmap=nothing)
    if ndims(A) < ndims(B)
        broadcast((x,y)->f(y,x), B, A; dimmap=dimmap)
    end

    if dimmap === nothing
        dimmap = match_indices(A, B)
    end

    common_A = Iterators.filter(i->dimmap[i] > 0, 1:ndims(A)) |> collect
    common_B = Iterators.filter(i -> i>0, dimmap) |> collect
    #@assert length(common_B) == length(common_A)

    # for every bounding box in A, take compare common_A
    # bounding boxes vs. every common_B bounding box

    out_chunks = []
    innerbcast(a, b) = broadcast(f, a, b; dimmap=dimmap)
    out_domains = []
    for (dA, cA) in zip(A.domains, A.chunks)
        for (dB, cB) in zip(B.domains, B.chunks)
            boxA = subbox(dA.boundingrect, common_A)
            boxB = subbox(dB.boundingrect, common_B)

            fst = map(max, boxA.first, boxB.first)
            lst = map(min, boxA.last, boxB.last)
            dmn = bcast_narrow_space(dA, common_A, fst, lst)
            if boxhasoverlap(boxA, boxB)
                push!(out_chunks, delayed(innerbcast)(cA, cB))
                push!(out_domains, dmn)
            end
        end
    end
    with_overlaps(out_domains, out_chunks) do chunks
        treereduce(delayed(_merge), chunks)
    end
end

broadcasted(f, A::DNDSparse, B::DNDSparse) = broadcast(f, A, B)

function match_indices(A::DNDSparse{K1},B::DNDSparse{K2}) where {K1,K2}
    if K1 <: NamedTuple && K2 <: NamedTuple
        Ap = dimlabels(A)
        Bp = dimlabels(B)
    else
        Ap = K1.parameters
        Bp = K2.parameters
    end
    IndexedTables.find_corresponding(Ap, Bp)
end

## Deprecation

Base.@deprecate naturaljoin(left::DNDSparse, right::DNDSparse, op::Function) naturaljoin(op, left::DNDSparse, right::DNDSparse)

Base.@deprecate leftjoin(left::DNDSparse, right::DNDSparse, op::Function) leftjoin(op, left, right)


text:
# What modules / packages do we depend on
using Random
using LinearAlgebra
using Printf
using Plots

# Import custom types and routines
using MyGEMM.Types: TilingParams
using MyGEMM.TransposeNoPacking: multiply_transpose_no_packing!
using MyGEMM.TransposeMultithreading: multiply_transpose!

Random.seed!(777) # To ensure repeatability
BLAS.set_num_threads(6)

# What precision numbers to use
FT = Float64

# Block matrix-matrix multiply parameters
if FT == Float64
    _mc = 96
    _nc = 2048
    _kc = 256
    _mr = 8
    _nr = 6
    _vl = 4
elseif FT == Float32
    _mc = 2 * 96
    _nc = 2 * 2048
    _kc = 2 * 256
    _mr = 16
    _nr = 12
    _vl = 8
end
params = TilingParams{_mc, _nc, _kc, _mr, _nr, _vl}()

# Different implementations
refgemm!(C, A, B) = mul!(C, A, B', one(eltype(C)), one(eltype(C)))
mygemm_no_packing!(C, A, B) = multiply_transpose_no_packing!(C, A, B, params::TilingParams)
mygemm_packing!(C, A, B) = multiply_transpose!(C, A, B, params::TilingParams)

# Benchmarking routine
function evaluate(gemm!, bench_dims, num_reps)
  num_runs = length(bench_dims)
  performances = []

  # Do loop over varying matrix sizes
  for nmk in bench_dims
    n = m = k = nmk
    gflops = 2 * m * n * k * 1e-09

    # Create some random initial data
    A = rand(FT, m, k)
    B = rand(FT, k, n)
    C = rand(FT, m, n)

    # Make a copy of C for resetting data later
    C_old = copy(C)

    # "truth"
    C_ref = C + A * B'

    # Compute the timings
    best_time = typemax(FT)
    for iter = 1:num_reps
      # Reset C to the original data
      C .= C_old;
      run_time = @elapsed gemm!(C, A, B);
      best_time = min(run_time, best_time)
    end

    # Make sure that we have the right answer!
    @assert C ≈ C_ref

    # Store data
    best_perf = gflops / best_time
    push!(performances, best_perf)
  end

  return performances
end

# Do performance evaluation
num_reps = 3
bench_dims = 48:48:2000
performances = []
for gemm! in [refgemm!, mygemm_no_packing!, mygemm_packing!]
  push!(performances, evaluate(gemm!, bench_dims, num_reps))
end

# Make simple performance plot
p0 = plot(bench_dims, performances[1], linewidth=2, label="reference", legend=:bottom, xlim=(minimum(bench_dims), maximum(bench_dims)))
plot!(p0, bench_dims, performances[2], linewidth=2, label="toby's GEMM no packing (nthreads=1)")
plot!(p0, bench_dims, performances[3], linewidth=2, label="toby's GEMM packing")
xlabel!(p0, "Matrix size")
ylabel!(p0, "GFLOPS")
title!(p0, "Performance evaluation (6 threads)")
savefig("performance_mt.png")

max_stars_repo_name: bischtob/MyGEMM

max_issues_repo_name: bischtob/MyGEMM

max_forks_repo_name: bischtob/MyGEMM


text:
"""
    lcs_length(str_1::String, str_2::String)

The longest common subsequence (LCS) algorithm is dedicated to solve the problem of finding
the longest subsequence common to all sequences in a set of two sequences. For finding these
subsequences are brute force approach (permuation) has to be chosen as follow:

    1. Figure out the number of possible different subsequences of a string with length `m` and `n`.
    2. Recontracted `for-loops` for defining the outer layer of the storing matrix `C`.
    3. Checking if `str_1[i] == str_2[j]` and updating the count in  matrix `C`.
    4.Return the latest element of `C[m,n]`.

The **LCS** differs from the [longest common substring problem](https://en.wikipedia.org/wiki/Longest_common_substring_problem) because it is not required to occupy consecutive positions within the original sequences.
For more information see: [https://en.wikipedia.org/wiki/Longest_common_subsequence_problem](https://en.wikipedia.org/wiki/Longest_common_subsequence_problem)

# Arguments
- `str_1::String`: String 1 of length `m`
- `str_2::String`: String 2 of length `n`

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> ClassicAlgorithmsCollections.lcs_lenght("AGGTAB","GXTXAYB")
4
```
"""
function lcs_length(str_1::String, str_2::String)
    m = length(str_1)
    n = length(str_2)
    C = zeros(Int64, (m, n))
    for i in 1:m
        C[i, 1] = 1
    end
    for j in 1:n
        C[1, j] = 1
    end
    for i in 2:m
        for j in 2:n
            if str_1[i] == str_2[j]
                C[i, j] = C[i-1, j-1] + 1
            else
                if C[i, j-1] > C[i-1, j]
                    C[i, j] = C[i, j-1]
                else
                    C[i, j] = C[i-1, j]
                end
            end
        end
    end
    return C[m, n]
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
# Vertex interface
"""
    vertex_id(vertex)

Return an identifier of the vertex. The identifier must be of a type that can be converted from and to an Int.
"""
function vertex_id end

"""
    set_vertex_id!(vertex, id)

Set the vertex identifier. The identifier must be of a type that can be converted from and to an Int.
"""
function set_vertex_id! end

"""
    vertex_id_type(V)

Return the identifier type used by vertex type `V`. The identifier type must be convertible from and to an Int.
"""
vertex_id_type(::Type) = Int # fallback, can be specialized

@inline vertex_index(x) = Int(vertex_id(x))
@inline vertex_index!(x::T, index::Int) where {T} = set_vertex_id!(x, vertex_id_type(T)(index))

# Edge interface
"""
    edge_id(edge)

Return an identifier of the edge. The identifier must be of a type that can be converted from and to an Int.
"""
function edge_id end

"""
    set_edge_id!(edge, id)

Set the edge identifier. The identifier must be of a type that can be converted from and to an Int.
"""
function set_edge_id! end

"""
    edge_id_type(V)

Return the identifier type used by edge type `V`. The identifier type must be convertible from and to an Int.
"""
edge_id_type(::Type) = Int # fallback, can be specialized

@inline edge_index(x) = Int(edge_id(x))
@inline edge_index!(x::T, index::Int) where {T} = set_edge_id!(x, edge_id_type(T)(index))

flip_direction(x::Any) = deepcopy(x)

abstract type AbstractGraph{V, E} end
vertextype(::Type{<:AbstractGraph{V, E}}) where {V, E} = V
vertextype(g::AbstractGraph) = vertextype(typeof(g))
edgetype(::Type{<:AbstractGraph{V, E}}) where {V, E} = E
edgetype(g::AbstractGraph) = edgetype(typeof(g))
num_vertices(g::AbstractGraph) = length(vertices(g))
num_edges(g::AbstractGraph) = length(edges(g))
out_neighbors(vertex::V, g::AbstractGraph{V, E}) where {V, E} = (target(e, g) for e in out_edges(vertex, g))
in_neighbors(vertex::V, g::AbstractGraph{V, E}) where {V, E} = (source(e, g) for e in in_edges(vertex, g))


text:
using JuliaTutorial
using Test

@testset "JuliaTutorial.jl" begin
    # @test myaddition(1, 1) == 2
end


text:
function recompile()
    curpath = pwd();
    cd(CONT_LIB_SRC)
    run(`make clean`)
    run(`make run CFLAGS='-fPIC -shared -DTRACE_LEVEL=2'`)
    cd(curpath)
end


text:
module PackageTest

include("extra_file.jl")
greet()  = println("Hello world")
println("Hello world")

end


text:
#
# Helper Functions
#

function readdata(filepath::AbstractString, sheet::Union{AbstractString, Int}, ref)
    xf = openxlsx(filepath, enable_cache=false)
    c = getdata(getsheet(xf, sheet), ref )
    close(xf)
    return c
end

function readdata(filepath::AbstractString, sheetref::AbstractString)
    xf = openxlsx(filepath, enable_cache=false)
    c = getdata(xf, sheetref)
    close(xf)
    return c
end

"""
    readtable(filepath, sheet, [columns]; [first_row], [column_labels], [header], [infer_eltypes], [stop_in_empty_row], [stop_in_row_function]) -> data, column_labels

Returns tabular data from a spreadsheet as a tuple `(data, column_labels)`.
`data` is a vector of columns. `column_labels` is a vector of symbols.
Use this function to create a `DataFrame` from package `DataFrames.jl`.

Use `columns` argument to specify which columns to get.
For example, `columns="B:D"` will select columns `B`, `C` and `D`.
If `columns` is not given, the algorithm will find the first sequence
of consecutive non-empty cells.

Use `first_row` to indicate the first row from the table.
`first_row=5` will look for a table starting at sheet row `5`.
If `first_row` is not given, the algorithm will look for the first
non-empty row in the spreadsheet.

`header` is a `Bool` indicating if the first row is a header.
If `header=true` and `column_labels` is not specified, the column labels
for the table will be read from the first row of the table.
If `header=false` and `column_labels` is not specified, the algorithm
will generate column labels. The default value is `header=true`.

Use `column_labels` as a vector of symbols to specify names for the header of the table.

Use `infer_eltypes=true` to get `data` as a `Vector{Any}` of typed vectors.
The default value is `infer_eltypes=false`.

`stop_in_empty_row` is a boolean indicating wether an empty row marks the end of the table.
If `stop_in_empty_row=false`, the `TableRowIterator` will continue to fetch rows until there's no more rows in the Worksheet.
The default behavior is `stop_in_empty_row=true`.

`stop_in_row_function` is a Function that receives a `TableRow` and returns a `Bool` indicating if the end of the table was reached.

Example for `stop_in_row_function`:

```
function stop_function(r)
    v = r[:col_label]
    return !Missings.ismissing(v) && v == "unwanted value"
end
```

Rows where all column values are equal to `Missing.missing` are dropped.

Example code for `readtable`:

```julia
julia> using DataFrames, XLSX

julia> df = DataFrame(XLSX.readtable("myfile.xlsx", "mysheet")...)

See also: `gettable`.
```
"""
function readtable(filepath::AbstractString, sheet::Union{AbstractString, Int}; first_row::Int = 1, column_labels::Vector{Symbol}=Vector{Symbol}(), header::Bool=true, infer_eltypes::Bool=false, stop_in_empty_row::Bool=true, stop_in_row_function::Union{Function, Void}=nothing, enable_cache::Bool=false)
    xf = openxlsx(filepath, enable_cache=enable_cache)
    c = gettable(getsheet(xf, sheet); first_row=first_row, column_labels=column_labels, header=header, infer_eltypes=infer_eltypes, stop_in_empty_row=stop_in_empty_row, stop_in_row_function=stop_in_row_function)
    close(xf)
    return c
end

"""
    writetable(filename, data, columnnames; [rewrite], [sheetname])

`data` is a vector of columns.
`columnames` is a vector of column labels.
`rewrite` is a `Bool` to control if `filename` should be rewrited if already exists.
`sheetname` is the name for the worksheet.

Example using `DataFrames.jl`:

```julia
import DataFrames, XLSX
df = DataFrames.DataFrame(integers=[1, 2, 3, 4], strings=["Hey", "You", "Out", "There"], floats=[10.2, 20.3, 30.4, 40.5])
XLSX.writetable("df.xlsx", DataFrames.columns(df), DataFrames.names(df))
```
"""
function writetable(filename::AbstractString, data, columnnames; rewrite::Bool=false, sheetname::AbstractString="", anchor_cell::Union{String, CellRef}=CellRef("A1"))

    if !rewrite
        @assert !isfile(filename) "$filename already exists."
    end

    xf = open_default_template()
    sheet = xf[1]

    if sheetname != ""
        rename!(sheet, sheetname)
    end

    if isa(anchor_cell, String)
        anchor_cell = CellRef(anchor_cell)
    end

    writetable!(sheet, data, columnnames; anchor_cell=anchor_cell)

    # write output file
    writexlsx(filename, xf, rewrite=rewrite)
    nothing
end

"""
    writetable(filename::AbstractString; rewrite::Bool=false, kw...)

Write multiple tables.

`kw` is a variable keyword argument list. Each element should be in this format: `sheetname=( data, column_names )`,
where `data` is a vector of columns and `column_names` is a vector of column labels.

Example:

```julia
import DataFrames, XLSX

df1 = DataFrames.DataFrame(COL1=[10,20,30], COL2=["Fist", "Sec", "Third"])
df2 = DataFrames.DataFrame(AA=["aa", "bb"], AB=[10.1, 10.2])

XLSX.writetable("report.xlsx", REPORT_A=( DataFrames.columns(df1), DataFrames.names(df1) ), REPORT_B=( DataFrames.columns(df2), DataFrames.names(df2) ))
```
"""
function writetable(filename::AbstractString; rewrite::Bool=false, kw...)
    if !rewrite
        @assert !isfile(filename) "$filename already exists."
    end

    xf = open_default_template()

    is_first = true

    for (sheetname, (data, column_names)) in kw
        if is_first
            # first sheet already exists in template file
            sheet = xf[1]
            rename!(sheet, string(sheetname))
            writetable!(sheet, data, column_names)

            is_first = false
        else
            sheet = addsheet!(xf, string(sheetname))
            writetable!(sheet, data, column_names)
        end
    end

    # write output file
    writexlsx(filename, xf, rewrite=rewrite)
    nothing
end


text:
using Documenter, ProximalOperators

makedocs(
  modules = [ProximalOperators],
  format = :html,
  sitename = "ProximalOperators.jl",
  authors = "Lorenzo Stella and contributors",
  pages = Any[
  "Home" => "index.md",
  "Functions" => "functions.md",
  "Calculus rules" => "calculus.md",
  "Prox and gradient" => "operators.md",
  "Demos" => "demos.md"
  ],
)

deploydocs(
  repo   = "github.com/kul-forbes/ProximalOperators.jl.git",
  julia  = "0.6",
  osname = "linux",
  target = "build",
  deps = nothing,
  make = nothing,
)


text:
export absolute_error

function Base.:(|>)(t::Transform{T}, p::Point{3,U}) where {T,U}
    (xₚ, yₚ, zₚ, wₚ) = sum(transpose(t.m[:,StaticArrays.SUnitRange(1,3)]) .* p,), dims=1) .+ transpose(t.m[:,StaticArrays.SUnitRange(4,4))]
    return wₚ ≈ 1 ? promote_type(T, typeof(p))(xₚ, yₚ, zₚ) / wₚ : promote_type(T, typeof(p))(xₚ, yₚ, zₚ)
end

absolute_error(t::Transform{T}, p::Point3{U}) where {T,U} =
    Γ(3.0) * Vector3{promote_type(T,U)}(
        sum(abs.(transpose(t.m[StaticArrays.SUnitRange(1,3), StaticArrays.SUnitRange(1,3)]) .* p), dims=1))

Base.:(|>)(t::Transform{T}, v::Vector3{U}) where {T,U} =
    Vector3{promote_type(T,U)}(
        sum(transpose(t.m[StaticArrays.SUnitRange(1,3), StaticArrays.SUnitRange(1,3)]) .* v, dims=1))

Base.:(|>)(t::Transform{T}, n::Normal3{U}}) where {T<:AbstractFloat, U<:Number} =
    Normal3{promote_type(T,U)}(
        sum(transpose(t.m_inv[StaticArrays.SUnitRange(1,3), StaticArrays.SUnitRange(1,3)]) .* n, dims=1))

function Base.:(|>)(t::Transform{T}, r::Ray{U}) where {T<:AbstractFloat, U<:AbstractFloat}
    m² = magnitude²(r.d)
    if m² > 0.0
        δₜ = (abs(r.d) ⋅ absolute_error(t, r.o)) / m²
        o = r.o + r.d * δₜ
        return Ray{promote_type(T,U)}(t |> o, t |> r.d, r.time)
    end
    Ray{promote_type(T,U)}(t |> r.o, t |> r.d, r.time)
end

function Base.:(|>)(t::Transform{T}, b::Bounds{N,U}) where {N,T,U}
    Bounds{N, promote_type(T,U)}(
        sum(min.(
            t.m[StaticArrays.SUnitRange(1,3), StaticArrays.SUnitRange(1,4)] .* b.p_min,
            t.m[StaticArrays.SUnitRange(1,3), StaticArrays.SUnitRange(1,4)] .* b.p_max)),
        sum(max.(
            t.m[StaticArrays.SUnitRange(1,3), StaticArrays.SUnitRange(1,4)] .* b.p_min,
            t.m[StaticArrays.SUnitRange(1,3), StaticArrays.SUnitRange(1,4)] .* b.p_max)))
end


text:
export AlternateSelector

"""
    AlternateSelector(n::Int)

Used to ensure that all actions are selected alternatively.

# Fields

- `n::Int`: means the optional actions are `1:n`.
- `step::Int=0`: record the number of times that the selector is applied.
"""
Base.@kwdef mutable struct AlternateSelector <: AbstractDiscreteActionSelector
    n::Int
    step::Int = 0
end

"""
    (s::AlternateSelector)(values::Any)

Here the `values` is ignored. The returned action is based on the `step` of `s`.

# Example

```julia
julia> selector = AlternateSelector(n=3)
AlternateSelector(3, 0)

julia> any_state = 0 # for AlternateSelector, the value of actions can be anything

julia> [selector(any_state) for i in 1:10]  # iterate through all actions
10-element Array{Int64,1}:
 1
 2
 3
 1
 ⋮
 1
 2
 3
 1
```
"""
function (s::AlternateSelector)(values::Any; kw...)
    s.step += 1
    (s.step - 1) % s.n + 1
end


text:
# Finding X-points and O-points in a 2D magnetic reconnection configuration.
# This example assumes X-Z meridional plane.
# Note:
# 1. The input B field domain matters for computing the flux function, but I'm entirely sure
# why there are differences.
# 2. In identifying the X-points and O-points, we currently provide two methods: method 1
# needs to set the relative tolerance, while method 2 does not.
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, PyPlot

file = "bulk.0001657.vlsv"

meta = load(file)

ndims(meta) != 2 && @error "Flux function only works for 2D simulations!"

if meta.ncells[3] == 1
   @error "equatorial plane, no reconnection!"
end

b = meta["vg_b_vol"]
b = reshape(b, 3, meta.ncells[1], meta.ncells[3])
dx = [meta.dcoord[1], meta.dcoord[3]]

x = LinRange(meta.coordmin[1], meta.coordmax[1], meta.ncells[1]) ./ Vlasiator.RE
z = LinRange(meta.coordmin[3], meta.coordmax[3], meta.ncells[3]) ./ Vlasiator.RE

xmin_ = searchsortedfirst(x, 7.0)
xmax_ = searchsortedlast(x, 9.5)
zmin_ = searchsortedfirst(z, -5.0)
zmax_ = searchsortedlast(z, 5.0)

# meshgrid for plotting
X = [a for a in x[xmin_:xmax_], _ in z[zmin_:zmax_]]
Z = [b for _ in x[xmin_:xmax_], b in z[zmin_:zmax_]]

flux = compute_flux_function(view(b,:,xmin_:xmax_,zmin_:zmax_), dx)

fig, ax = plt.subplots(subplot_kw=Dict("projection"=>"3d"))

ax.plot_surface(X, Z, flux;
   cmap=matplotlib.cm.turbo,
   linewidth=0, antialiased=false)

indices_x, indices_o = find_reconnection_points(flux; retol=1e-4, method=2)

fig, ax = subplots(figsize=(6,10), constrained_layout=true)

pcolormesh(meta, "proton/vg_v", ax;
   comp=:z, extent=[5, 10, -7.5, 7.5],
   cmap=matplotlib.cm.RdBu_r)
s1 = ax.scatter(x[indices_x[1,:].+xmin_.-1], z[indices_x[2,:].+zmin_.-1];
   s=50, marker="x", color="tab:gray")
s2 = ax.scatter(x[indices_o[1,:].+xmin_.-1], z[indices_o[2,:].+zmin_.-1];
   s=12, marker="o", color="tab:brown")

ax.legend([s1, s2], ["X-point", "O-point"])


text:
export find_minimum, roots

# Quadratic Polynomial to represent cost-to-go function,
# keeos track of its "ancestor" to simplify the recovery of the solution
# has_been_used keeps track of if the Quadratic polynimial has been inserted
# into a piecewise quadratic, otherwise it is reused
mutable struct QuadraticPolynomial{T<:Real}
a::T
b::T
c::T
has_been_used::Bool
time_index::Int32
ancestor::QuadraticPolynomial{T}
function QuadraticPolynomial{T}(a::T, b::T, c::T) where {T}
    # @assert a ≥ 0, # Δ may have negative a ..., Δ is currently not used...
    new(a, b, c, false, -1)
end

QuadraticPolynomial{T}() where {T} = new()
end

QuadraticPolynomial(a::T,b::T,c::T) where T = QuadraticPolynomial{T}(a,b,c)
QuadraticPolynomial(a,b,c) = QuadraticPolynomial(promote(a,b,c)...)

function QuadraticPolynomial(x::Vector{T}) where T
    QuadraticPolynomial(x[3], x[2], x[1])
end

function Base.show(io::IO, p::QuadraticPolynomial)
    @printf(io, "%.2f*x^2 + %.2f*x + %.2f   ", p.a, p.b, p.c)
end

function Base.zero(::Type{QuadraticPolynomial{T}}) where T
    return QuadraticPolynomial(zero(T), zero(T), zero(T))
end
Base.zero(S::T) where {T<:QuadraticPolynomial} = zero(T)

-(p1::QuadraticPolynomial, p2::QuadraticPolynomial) = QuadraticPolynomial(p1.a-p2.a, p1.b-p2.b, p1.c-p2.c)

==(p1::QuadraticPolynomial, p2::QuadraticPolynomial) = (p1.a==p2.a) && (p1.b==p2.b) && (p1.c==p2.c)

@inline function (p::QuadraticPolynomial)(x)
    return p.a*x^2 + p.b*x + p.c
end

# Finds the minimum of a positive definite quadratic one variable polynomial
# the find_minimum fcn returns (opt_x, opt_val)
function find_minimum(p::QuadraticPolynomial)
    if p.a < 0 || (p.a == 0 && p.b != 0)
        println("No unique minimum exists")
        return (NaN, NaN)
    elseif p.a == 0
        # p.b == 0
        return (NaN, p.c)
    else
        x_opt = -p.b / 2 / p.a
        f_opt = -p.b^2/4/p.a + p.c
        return (x_opt, f_opt)
    end
end

"""
Minimum of a quadratic function which is assumed to be positive definite
No checks of this is done
"""
@inline function unsafe_minimum(p::QuadraticPolynomial{T}) where T
    return (-p.b^2/4/p.a + p.c)::T
end


text:
using Revise, Lifting
import Lifting:
    Lifting_Programmes, Lifting_Aux, Lifting_Exercise_Names, Lifting_Progressions
cd(@__DIR__)
name = "nSunsCAP3_OHP_6Day_LP"

# Uncomment to print whole programme.
println(Lifting_Programmes[name])
write("programme.csv", Lifting_Programmes[name])
Lifting_Programmes[name].exerProg["bulgarianSplit"].progression
Lifting_Programmes[name]

calcIntensity(6, 8.5)

println(keys(Lifting_Programmes[name].exerProg))

name = "DCG_YugeStr"
Lifting_Programmes[name].days
println(Lifting_Programmes[name], 1)

write("Yuge_n_Stronk.csv", Lifting_Programmes[name])

max_stars_repo_name: dcelisgarza/Lifting


text:
using Test
using ParticleInCell

@testset "ParticleMeshCoupling 2D" begin

    n_cells = 10 # Number of cells
    n_particles = 4 # Number of particles
    spline_degree = 3 # Spline degree
    xmin, xmax, nx = 0.0, 2.0, n_cells
    ymin, ymax, ny = 0.0, 1.0, n_cells

    mesh = TwoDGrid(xmin, xmax, nx, ymin, ymax, ny)

    volume = (xmax - xmin) * (ymax - ymin)

    x_vec = [0.1 0.65 0.7 1.5; 0.0 0.0 0.0 0.0]
    v_vec = [1.5 0.00 0.0 0.0; 0.0 0.5 0.0 0.0]

    pg = ParticleGroup{2,2}(n_particles, charge = 1.0, mass = 1.0, n_weights = 1)

    for i_part = 1:n_particles
        pg.array[1, i_part] = x_vec[1, i_part]
        pg.array[2, i_part] = x_vec[2, i_part]
        pg.array[3, i_part] = v_vec[1, i_part]
        pg.array[4, i_part] = v_vec[2, i_part]
        pg.array[5, i_part] = 1.0 / n_particles
    end

    # Initialize the kernel
    kernel = ParticleMeshCoupling2D(pg, mesh, spline_degree, :collocation)

    # Reference values of the shape factors
    index_grid = zeros(Int, (2, 4))
    index_grid[1, :] .= [-2, 1, 1, 5]
    index_grid[2, :] .= [-3, -3, -3, -3]

    v_grid = zeros(Float64, (4, 2, 4))
    v_grid[:, 1, 1] .= [
        2.0833333333333332E-002,
        0.47916666666666663,
        0.47916666666666663,
        2.0833333333333332E-002,
    ]
    v_grid[:, 1, 3] .= v_grid[:, 1, 1]
    v_grid[:, 1, 4] .= v_grid[:, 1, 1]
    v_grid[:, 1, 2] .= [
        7.0312500000000000E-002,
        0.61197916666666663,
        0.31510416666666663,
        2.6041666666666665E-003,
    ]
    v_grid[1, 2, :] .= 0.0
    v_grid[2, 2, :] .= 1.0 / 6.0
    v_grid[3, 2, :] .= 2.0 / 3.0
    v_grid[4, 2, :] .= 1.0 / 6.0

    ρ_dofs = zeros(Float64, 100)
    ρ_dofs1 = zeros(Float64, 100)
    ρ_dofs_ref = zeros(Float64, 100)
    ρ_dofs_pp = zeros(Float64, (16, 100))

    # Accumulate ρ
    for i_part = 1:n_particles
        x = pg.array[1, i_part]
        y = pg.array[2, i_part]
        w = pg.array[5, i_part]
        GEMPIC.add_charge!(ρ_dofs, kernel, x, y, w)
        GEMPIC.add_charge_pp!(ρ_dofs1, kernel, x, y, w)
    end

    ρ_dofs_ref[8:10] .= v_grid[1:3, 1, 1]
    ρ_dofs_ref[1] = v_grid[4, 1, 1]
    ρ_dofs_ref[1:4] .= ρ_dofs_ref[1:4] .+ v_grid[:, 1, 2] .+ v_grid[:, 1, 3]
    ρ_dofs_ref[5:8] .= ρ_dofs_ref[5:8] .+ v_grid[:, 1, 4]

    ρ_dofs_ref[71:80] .= ρ_dofs_ref[1:10] ./ 6.0
    ρ_dofs_ref[81:90] .= ρ_dofs_ref[1:10] .* 2.0 / 3.0
    ρ_dofs_ref[91:100] .= ρ_dofs_ref[1:10] ./ 6.0
    ρ_dofs_ref[1:10] .= 0.0

    ρ_dofs_ref .*= n_cells^2 / volume / n_particles

    @test ρ_dofs ≈ ρ_dofs_ref
    @test ρ_dofs1 ≈ ρ_dofs_ref

    GEMPIC.b_to_pp_2d!(ρ_dofs_pp, kernel.spline1, kernel.spline2, ρ_dofs)

    # For evaluation check
    particle_values = zeros(4)
    particle_values1 = zeros(4)
    particle_values_ref = zeros(4)

    # Test function evaluation
    for i_part = 1:n_particles
        xp = pg.array[1, i_part]
        yp = pg.array[2, i_part]
        particle_values[i_part] = evaluate(kernel, xp, yp, ρ_dofs)
        particle_values1[i_part] = evaluate_pp(kernel, xp, yp, ρ_dofs_pp)
    end

    particle_values_ref =
        [1.1560058593749998, 2.3149278428819446, 2.2656250000000000, 1.1512586805555554] ./
        volume

    fill!(particle_values_ref, 0.0)

    for i_part = 1:n_particles
        for i = 1:4
            i1 = mod(index_grid[1, i_part] + i - 2, n_cells)
            for j = 1:4
                i2 = mod(index_grid[2, i_part] + j - 2, n_cells)
                res =
                    v_grid[i, 1, i_part] *
                    v_grid[j, 2, i_part] *
                    ρ_dofs_ref[i1+i2*n_cells+1]
                particle_values_ref[i_part] += res
            end
        end
    end

    @test particle_values ≈ particle_values_ref
    @test particle_values1 ≈ particle_values_ref

end


text:
using LowRankModels
using FactCheck

TOL = 1e-3

facts("Classification Losses") do

  context("logistic loss") do
    l = LogisticLoss()
    l1 = 1.31326168
    l0 = 0.3132616875
    # 1 is true
    # -1 and 0 are both false
    # anything else is an error
    @fact evaluate(l, 1, true) --> roughly(l0, TOL)
    @fact evaluate(l, 1, false) --> roughly(l1, TOL)
    @fact evaluate(l, -1, true) --> roughly(l1, TOL)
    @fact evaluate(l, -1, false) --> roughly(l0, TOL)
    @fact evaluate(l, 1, 1) --> roughly(l0, TOL)
    @fact evaluate(l, 1, -1) --> roughly(l1, TOL)
    @fact evaluate(l, 1, 0) --> roughly(l1, TOL)
    @fact evaluate(l, -1, 1) --> roughly(l1, TOL)
    @fact evaluate(l, -1, -1) --> roughly(l0, TOL)
    @fact evaluate(l, -1, 0) --> roughly(l0, TOL)

    @fact evaluate(3*l, 1, false) --> roughly(3*l1, TOL)
  end

  context("hinge loss") do
    l = HingeLoss()
    # 1 is true
    # -1 and 0 are both false
    # anything else is an error
    @fact evaluate(l, 1, true) --> roughly(0, TOL)
    @fact evaluate(l, 1, false) --> roughly(2, TOL)
    @fact evaluate(l, -1, true) --> roughly(2, TOL)
    @fact evaluate(l, -1, false) --> roughly(0, TOL)
    @fact evaluate(l, 1, 1) --> roughly(0, TOL)
    @fact evaluate(l, 1, -1) --> roughly(2, TOL)
    @fact evaluate(l, 1, 0) --> roughly(2, TOL)
    @fact evaluate(l, -1, 1) --> roughly(2, TOL)
    @fact evaluate(l, -1, -1) --> roughly(0, TOL)
    @fact evaluate(l, -1, 0) --> roughly(0, TOL)

    @fact evaluate(3*l, 1, false) --> roughly(3*2, TOL)

    @fact grad(l, -1, true) --> roughly(-1, TOL)
    @fact grad(l, 2, true) --> roughly(0, TOL)
    @fact grad(l, -2, false) --> roughly(0, TOL)
    @fact grad(l, 2, false) --> roughly(1, TOL)
  end

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 ApplicationGatewayBackendHealth <: SwaggerModel
    backendAddressPools::Any # spec type: Union{ Nothing, Vector{ApplicationGatewayBackendHealthPool} } # spec name: backendAddressPools

    function ApplicationGatewayBackendHealth(;backendAddressPools=nothing)
        o = new()
        validate_property(ApplicationGatewayBackendHealth, Symbol("backendAddressPools"), backendAddressPools)
        setfield!(o, Symbol("backendAddressPools"), backendAddressPools)
        o
    end
end # type ApplicationGatewayBackendHealth

const _property_map_ApplicationGatewayBackendHealth = Dict{Symbol,Symbol}(Symbol("backendAddressPools")=>Symbol("backendAddressPools"))
const _property_types_ApplicationGatewayBackendHealth = Dict{Symbol,String}(Symbol("backendAddressPools")=>"Vector{ApplicationGatewayBackendHealthPool}")
Base.propertynames(::Type{ ApplicationGatewayBackendHealth }) = collect(keys(_property_map_ApplicationGatewayBackendHealth))
Swagger.property_type(::Type{ ApplicationGatewayBackendHealth }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayBackendHealth[name]))}
Swagger.field_name(::Type{ ApplicationGatewayBackendHealth }, property_name::Symbol) =  _property_map_ApplicationGatewayBackendHealth[property_name]

function check_required(o::ApplicationGatewayBackendHealth)
    true
end

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


text:
module GA
    using Polyester, LoopVectorization, VectorizedRNG

    function fit(;models::Array, input_data::Array{Float32}, output_data::Array{Float32}, monitor::Any, α::Float64=0.01, num_copy::Int64, epochs::Int64=20, batch::Real=32, mini_batch::Int64=5)
        gene_pool = length(models)
        input_shape = models[1].layers[2].input_shape
        output_shape = models[1].layers[end].output_shape
        examples = size(input_data)[end]
        batch_size = ceil(Int64, examples/batch)*mini_batch
        current_input_data = zeros(Float32, input_shape..., mini_batch)
        current_output_data = zeros(Float32, output_shape..., mini_batch)

        @batch for i in 1:gene_pool
            models[i].initialize(models[i], mini_batch)
        end

        for e in 1:epochs
            print("Epoch ", e, "\n[")
            loss = 0
            losses = zeros(Float32, gene_pool)
            weights = zeros(Float32, gene_pool)
            @time begin

                for t in 1:mini_batch:batch_size-mini_batch+1
                    @batch for i in 1:mini_batch
                        index = rand(1:examples)
                        selectdim(current_input_data, length(input_shape)+1, i) .= selectdim(input_data, length(input_shape)+1, index)
                        selectdim(current_output_data, length(output_shape)+1, i) .= selectdim(output_data, length(output_shape)+1, index)
                    end

                    if ((t+mini_batch-1)÷mini_batch)%ceil(batch_size/(50*mini_batch))==0
                        print("=")
                    end

                    for i in 1:gene_pool
                        models[i].activate(models[i], current_input_data)
                        losses[i] = monitor.func(models[i].layers[end-1].output, current_output_data)
                    end

                    for i in 1:gene_pool-num_copy
                        get_weights!(weights, losses)
                        recomutation!(models[argmax(losses)], models[sample(weights)], models[sample(weights)], α, t, batch_size-mini_batch+1)
                        losses[argmax(losses)] = Inf32
                    end

                    loss += minimum(losses)
                end
                print("] with loss ", loss, ", time usage ")
            end
        end
    end

    function get_weights!(weights::Array{Float32}, losses::Array{Float32})
        @avxt for i in eachindex(losses)
            weights[i] = ifelse(losses[i]!=Inf32, 1/losses[i], 0.0f0)
        end
        s = sum(weights)
        @avxt weights ./= s
    end

    function sample(weights)
        r = rand()
        for i in 1:length(weights)
            if weights[i]>=r
                return i
            else
                r -= weights[i]
            end
        end
    end

    function mutation_func(t, v, T)
        if rand(-1:2:1)>0
            return (1.0f0-v)*(1.0f0-rand(0.0f0:1.0f-3:1.0f0)^(1-t/T)^5)
        else
            return -(v+1.0f0)*(1.0f0-rand(0.0f0:1.0f-3:1.0f0)^(1-t/T)^5)
        end
    end

    function recomutation!(new_model, model₁, model₂, α, t, T)
        for i in 1:model₁.num_layer
            if hasproperty(new_model.layers[i], :filters)
                @avxt temp = model₂.layers[i].filters .- model₁.layers[i].filters
                rand!(local_rng(), new_model.layers[i].filters, VectorizedRNG.StaticInt(0), model₁.layers[i].filters, temp)
                rand!(local_rng(), temp)
                @avxt for j in eachindex(temp)
                    new_model.layers[i].filters[j] = ifelse(temp[j]<=α, ifelse(temp[j]<=0.5, (1.0f0-new_model.layers[i].filters[j])*(1.0f0-temp[j]^(1-t/T)^5), -(new_model.layers[i].filters[j]+1.0f0)*(1.0f0-temp[j]^(1-t/T)^5)), new_model.layers[i].filters[j])
                end
            end

            if hasproperty(new_model.layers[i], :weights)
                @avxt temp = model₂.layers[i].weights .- model₁.layers[i].weights
                rand!(local_rng(), new_model.layers[i].weights, VectorizedRNG.StaticInt(0), model₁.layers[i].weights, temp)
                rand!(local_rng(), temp)
                @avxt for j in eachindex(temp)
                    new_model.layers[i].weights[j] = ifelse(temp[j]<=α, ifelse(temp[j]<=0.5, (1.0f0-new_model.layers[i].weights[j])*(1.0f0-temp[j]^(1-t/T)^5), -(new_model.layers[i].weights[j]+1.0f0)*(1.0f0-temp[j]^(1-t/T)^5)), new_model.layers[i].weights[j])
                end
            end

            if hasproperty(new_model.layers[i], :biases)
                @avxt temp = model₂.layers[i].biases .- model₁.layers[i].biases
                rand!(local_rng(), new_model.layers[i].biases, VectorizedRNG.StaticInt(0), model₁.layers[i].biases, temp)
                rand!(local_rng(), temp)
                @avxt for j in eachindex(temp)
                    new_model.layers[i].biases[j] = ifelse(temp[j]<=α, ifelse(temp[j]<=0.5, (1.0f0-new_model.layers[i].biases[j])*(1.0f0-temp[j]^(1-t/T)^5), -(new_model.layers[i].biases[j]+1.0f0)*(1.0f0-temp[j]^(1-t/T)^5)), new_model.layers[i].biases[j])
                end
            end
        end
    end
end


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

using Test
isdefined(Main, :OffsetArrays) || @eval Main include("testhelpers/OffsetArrays.jl")
using .Main.OffsetArrays
isdefined(Main, :TSlow) || @eval Main include("testhelpers/arrayindexingtypes.jl")
using .Main: TSlow, WrapperArray

A = Int64[1, 2, 3, 4]
As = TSlow(A)
Ars = Int64[1 3; 2 4]
Arss = TSlow(Ars)
B = Complex{Int64}[5+6im, 7+8im, 9+10im]
Bs = TSlow(B)
Av = [Int32[1,2], Int32[3,4]]

for Ar in (Ars, Arss)
    @test @inferred(ndims(reinterpret(reshape, Complex{Int64}, Ar))) == 1
    @test @inferred(axes(reinterpret(reshape, Complex{Int64}, Ar))) === (Base.OneTo(2),)
    @test @inferred(size(reinterpret(reshape, Complex{Int64}, Ar))) == (2,)
end
for _B in (B, Bs)
    @test @inferred(ndims(reinterpret(reshape, Int64, _B))) == 2
    @test @inferred(axes(reinterpret(reshape, Int64, _B))) === (Base.OneTo(2), Base.OneTo(3))
    @test @inferred(size(reinterpret(reshape, Int64, _B))) == (2, 3)
    @test @inferred(ndims(reinterpret(reshape, Int128, _B))) == 1
    @test @inferred(axes(reinterpret(reshape, Int128, _B))) === (Base.OneTo(3),)
    @test @inferred(size(reinterpret(reshape, Int128, _B))) == (3,)
end

@test_throws ArgumentError("cannot reinterpret `Int64` as `Vector{Int64}`, type `Vector{Int64}` is not a bits type") reinterpret(Vector{Int64}, A)
@test_throws ArgumentError("cannot reinterpret `Vector{Int32}` as `Int32`, type `Vector{Int32}` is not a bits type") reinterpret(Int32, Av)
@test_throws ArgumentError("cannot reinterpret a zero-dimensional `Int64` array to `Int32` which is of a different size") reinterpret(Int32, reshape([Int64(0)]))
@test_throws ArgumentError("cannot reinterpret a zero-dimensional `Int32` array to `Int64` which is of a different size") reinterpret(Int64, reshape([Int32(0)]))
@test_throws ArgumentError("""cannot reinterpret an `$Int` array to `Tuple{$Int, $Int}` whose first dimension has size `5`.
                              The resulting array would have non-integral first dimension.
                              """) reinterpret(Tuple{Int,Int}, [1,2,3,4,5])

@test_throws ArgumentError("`reinterpret(reshape, Complex{Int64}, a)` where `eltype(a)` is Int64 requires that `axes(a, 1)` (got Base.OneTo(4)) be equal to 1:2 (from the ratio of element sizes)") reinterpret(reshape, Complex{Int64}, A)
@test_throws ArgumentError("`reinterpret(reshape, T, a)` requires that one of `sizeof(T)` (got 24) and `sizeof(eltype(a))` (got 16) be an integer multiple of the other") reinterpret(reshape, NTuple{3, Int64}, B)
@test_throws ArgumentError("cannot reinterpret `Int64` as `Vector{Int64}`, type `Vector{Int64}` is not a bits type") reinterpret(reshape, Vector{Int64}, Ars)
@test_throws ArgumentError("cannot reinterpret a zero-dimensional `UInt8` array to `UInt16` which is of a larger size") reinterpret(reshape, UInt16, reshape([0x01]))

# getindex
for _A in (A, As)
    @test reinterpret(Complex{Int64}, _A) == [1 + 2im, 3 + 4im]
    @test reinterpret(Float64, _A) == reinterpret.(Float64, A)
    @test reinterpret(reshape, Float64, _A) == reinterpret.(Float64, A)
end
for Ar in (Ars, Arss)
    @test reinterpret(reshape, Complex{Int64}, Ar) == [1 + 2im, 3 + 4im]
    @test reinterpret(reshape, Float64, Ar) == reinterpret.(Float64, Ars)
end

for _B in (B, Bs)
    @test reinterpret(NTuple{3, Int64}, _B) == [(5,6,7),(8,9,10)]
    @test reinterpret(reshape, Int64, _B) == [5 7 9; 6 8 10]
end

# setindex
for (_A, Ar, _B) in ((A, Ars, B), (As, Arss, Bs))
    let Ac = copy(_A), Arsc = copy(Ar), Bc = copy(_B)
        reinterpret(Complex{Int64}, Ac)[2] = -1 - 2im
        @test Ac == [1, 2, -1, -2]
        reinterpret(Complex{Int64}, Arsc)[2] = -1 - 2im
        @test Arsc == [1 -1; 2 -2]
        reinterpret(NTuple{3, Int64}, Bc)[2] = (4,5,6)
        @test Bc == Complex{Int64}[5+6im, 7+4im, 5+6im]
        reinterpret(NTuple{3, Int64}, Bc)[1] = (1,2,3)
        @test Bc == Complex{Int64}[1+2im, 3+4im, 5+6im]
        Bc = copy(_B)
        Brrs = reinterpret(reshape, Int64, Bc)
        Brrs[2, 3] = -5
        @test Bc == Complex{Int64}[5+6im, 7+8im, 9-5im]
        Brrs[last(eachindex(Brrs))] = 22
        @test Bc == Complex{Int64}[5+6im, 7+8im, 9+22im]

        A1 = reinterpret(Float64, _A)
        A2 = reinterpret(ComplexF64, _A)
        A1[1] = 1.0
        @test real(A2[1]) == 1.0
        A1 = reinterpret(reshape, Float64, _A)
        A1[1] = 2.5
        @test reinterpret(Float64, _A[1]) == 2.5
        A1rs = reinterpret(Float64, Ar)
        A2rs = reinterpret(ComplexF64, Ar)
        A1rs[1, 1] = 1.0
        @test real(A2rs[1]) == 1.0
        A1rs = reinterpret(reshape, Float64, Ar)
        A2rs = reinterpret(reshape, ComplexF64, Ar)
        A1rs[1, 1] = 2.5
        @test real(A2rs[1]) == 2.5
    end
end
A3 = collect(reshape(1:18, 2, 3, 3))
A3r = reinterpret(reshape, Complex{Int}, A3)
@test A3r[4] === A3r[1,2] === A3r[CartesianIndex(1, 2)] === 7+8im
A3r[2,3] = -8-15im
@test A3[1,2,3] == -8
@test A3[2,2,3] == -15
A3r[4] = 100+200im
@test A3[1,1,2] == 100
@test A3[2,1,2] == 200
A3r[CartesianIndex(1,2)] = 300+400im
@test A3[1,1,2] == 300
@test A3[2,1,2] == 400

# same-size reinterpret where one of the types is non-primitive
let a = NTuple{4,UInt8}[(0x01,0x02,0x03,0x04)]
    @test reinterpret(Float32, a)[1] == reinterpret(Float32, 0x04030201)
    reinterpret(Float32, a)[1] = 2.0
    @test reinterpret(Float32, a)[1] == 2.0
end
let a = NTuple{4,UInt8}[(0x01,0x02,0x03,0x04)]
    @test reinterpret(reshape, Float32, a)[1] == reinterpret(Float32, 0x04030201)
    reinterpret(reshape, Float32, a)[1] = 2.0
    @test reinterpret(reshape, Float32, a)[1] == 2.0
end

# Pass-through indexing
B = Complex{Int64}[5+6im, 7+8im, 9+10im]
Br = reinterpret(reshape, Int64, B)
W = WrapperArray(Br)
for (b, w) in zip(5:10, W)
    @test b == w
end
for (i, j) in zip(eachindex(W), 11:16)
    W[i] = j
end
@test B[1] === Complex{Int64}(11+12im)
@test B[2] === Complex{Int64}(13+14im)
@test B[3] === Complex{Int64}(15+16im)
z3 = (0x00, 0x00, 0x00)
Az = [z3 z3; z3 z3]
Azr = reinterpret(reshape, UInt8, Az)
W = WrapperArray(Azr)
copyto!(W, fill(0x01, 3, 2, 2))
@test all(isequal((0x01, 0x01, 0x01)), Az)

# ensure that reinterpret arrays aren't erroneously classified as strided
let A = reshape(1:20, 5, 4)
    V = view(A, :, :)
    R = reinterpret(Int32, V)
    R2 = reinterpret(Int32, A)
    @test !(R isa StridedArray)
    @test !(R2 isa StridedArray)
    @test R * ones(4, 5) == R2 * ones(4,5) == copy(R) * ones(4,5) == copy(R2) * ones(4,5)
end

# but ensure that strided views of strided reinterpret arrays are still strided
let A = collect(reshape(1:20, 5, 4))
    R = reinterpret(Int32, A)
    @test R isa StridedArray
    @test view(R, :, :) isa StridedArray
    @test reshape(R, :) isa StridedArray
end

# and ensure a reinterpret array containing a strided array can have strides computed
let A = view(reinterpret(Int16, collect(reshape(UnitRange{Int64}(1, 20), 5, 4))), :, 1:2)
    R = reinterpret(Int32, A)
    @test strides(R) == (1, 10)
    @test stride(R, 1) == 1
    @test stride(R, 2) == 10
end

@testset "strides" begin
    a = rand(10)
    b = view(a,2:2:10)
    A = rand(10,10)
    B = view(A, 2:2:10, 2:2:10)

    @test strides(a) == (1,)
    @test strides(b) == (2,)
    @test strides(A) == (1,10)
    @test strides(B) == (2,20)

    for M in (a, b, A, B)
        @inferred strides(M)
        strides_M = strides(M)

        for (i, _stride) in enumerate(collect(strides_M))
            @test _stride == stride(M, i)
        end
    end
end

# IndexStyle
let a = fill(1.0, 5, 3)
    r = reinterpret(Int64, a)
    @test @inferred(IndexStyle(r)) == IndexLinear()
    fill!(r, 2)
    @test all(a .=== reinterpret(Float64, [Int64(2)])[1])
    @test all(r .=== Int64(2))
    for badinds in (0, 16, (0,1), (1,0), (6,3), (5,4))
        @test_throws BoundsError r[badinds...]
        @test_throws BoundsError r[badinds...] = -2
    end
    for goodinds in (1, 15, (1,1), (5,3))
        r[goodinds...] = -2
        @test r[goodinds...] == -2
    end
    r = reinterpret(Int32, a)
    @test @inferred(IndexStyle(r)) == IndexLinear()
    fill!(r, 3)
    @test all(a .=== reinterpret(Float64, [(Int32(3), Int32(3))])[1])
    @test all(r .=== Int32(3))
    for badinds in (0, 31, (0,1), (1,0), (11,3), (10,4))
        @test_throws BoundsError r[badinds...]
        @test_throws BoundsError r[badinds...] = -3
    end
    for goodinds in (1, 30, (1,1), (10,3))
        r[goodinds...] = -3
        @test r[goodinds...] == -3
    end
    r = reinterpret(Int64, view(a, 1:2:5, :))
    @test @inferred(IndexStyle(r)) == IndexCartesian()
    fill!(r, 4)
    @test all(a[1:2:5,:] .=== reinterpret(Float64, [Int64(4)])[1])
    @test all(r .=== Int64(4))
    for badinds in (0, 10, (0,1), (1,0), (4,3), (3,4))
        @test_throws BoundsError r[badinds...]
        @test_throws BoundsError r[badinds...] = -4
    end
    for goodinds in (1, 9, (1,1), (3,3))
        r[goodinds...] = -4
        @test r[goodinds...] == -4
    end
    r = reinterpret(Int32, view(a, 1:2:5, :))
    @test @inferred(IndexStyle(r)) == IndexCartesian()
    fill!(r, 5)
    @test all(a[1:2:5,:] .=== reinterpret(Float64, [(Int32(5), Int32(5))])[1])
    @test all(r .=== Int32(5))
    for badinds in (0, 19, (0,1), (1,0), (7,3), (6,4))
        @test_throws BoundsError r[badinds...]
        @test_throws BoundsError r[badinds...] = -5
    end
    for goodinds in (1, 18, (1,1), (6,3))
        r[goodinds...] = -5
        @test r[goodinds...] == -5
    end
end

# Error on reinterprets that would expose padding
struct S1
    a::Int8
    b::Int64
end

struct S2
    a::Int16
    b::Int64
end

A1 = S1[S1(0, 0)]
A2 = S2[S2(0, 0)]
@test reinterpret(S1, A2)[1] == S1(0, 0)
@test_throws Base.PaddingError (reinterpret(S1, A2)[1] = S2(1, 2))
@test_throws Base.PaddingError reinterpret(S2, A1)[1]
reinterpret(S2, A1)[1] = S2(1, 2)
@test A1[1] == S1(1, 2)

# Unconventional axes
let a = [0.1 0.2; 0.3 0.4], at = reshape([(i,i+1) for i = 1:2:8], 2, 2)
    v = OffsetArray(a, (-1, 1))
    r = reinterpret(Int64, v)
    @test axes(r) === axes(v)
    @test r[0,2] === reinterpret(Int64, v[0,2])
    @test r[1,2] === reinterpret(Int64, v[1,2])
    @test r[0,3] === reinterpret(Int64, v[0,3])
    @test r[1,3] === reinterpret(Int64, v[1,3])
    @test_throws ArgumentError("cannot reinterpret a `Float64` array to `UInt32` when the first axis is OffsetArrays.IdOffsetRange(0:1). Try reshaping first.") reinterpret(UInt32, v)
    @test_throws ArgumentError("`reinterpret(reshape, Tuple{Float64, Float64}, a)` where `eltype(a)` is Float64 requires that `axes(a, 1)` (got OffsetArrays.IdOffsetRange(0:1)) be equal to 1:2 (from the ratio of element sizes)") reinterpret(reshape, Tuple{Float64,Float64}, v)
    v = OffsetArray(a, (0, 1))
    @test axes(reinterpret(reshape, Tuple{Float64,Float64}, v)) === (OffsetArrays.IdOffsetRange(Base.OneTo(2), 1),)
    r = reinterpret(UInt32, v)
    axsv = axes(v)
    @test axes(r) === (oftype(axsv[1], 1:4), axsv[2])
    for i = 1:2
        rval = reinterpret(Tuple{UInt32,UInt32}, [v[i,2]])[1]
        @test r[2i-1,2]   == rval[1]
        @test r[2i,2] == rval[2]
        rval = reinterpret(Tuple{UInt32,UInt32}, [v[i,3]])[1]
        @test r[2i-1,3]   == rval[1]
        @test r[2i,3] == rval[2]
    end
    r[4,2] = 7
    @test r[4,2] === UInt32(7)
    @test a[2,1] === reinterpret(Float64, [0x33333333, UInt32(7)])[1]
    offsetvt = (-2, 4)
    vt = OffsetArray(at, offsetvt)
    istr = string(Int)
    @test_throws ArgumentError("cannot reinterpret a `Tuple{$istr, $istr}` array to `$istr` when the first axis is OffsetArrays.IdOffsetRange(-1:0). Try reshaping first.") reinterpret(Int, vt)
    vt = reshape(vt, 1:1, axes(vt)...)
    r = reinterpret(Int, vt)
    @test r == OffsetArray(reshape(1:8, 2, 2, 2), (0, offsetvt...))
end

@testset "potentially aliased copies" begin
    buffer = UInt8[1,0,0,0,0,0,0,0,2,0,0,0,0,0,0,0]
    mid = length(buffer) ÷ 2
    x1 = reinterpret(Int64, @view buffer[1:mid])
    x2 = reinterpret(Int64, @view buffer[mid+1:end])
    x1 .= x2
    @test x1 == x2 == [2]
    @test x1[] === x2[] === Int64(2)
end

# avoid nesting
@test parent(reinterpret(eltype(A), reinterpret(eltype(B), A))) === A

# Test 0-dimensional Arrays
A = zeros(UInt32)
B = reinterpret(Int32,A)
Brs = reinterpret(reshape,Int32,A)
@test size(B) == size(Brs) == ()
@test axes(B) == axes(Brs) == ()
B[] = Int32(5)
@test B[] === Int32(5)
@test Brs[] === Int32(5)
@test A[] === UInt32(5)

a = [(1.0,2.0)]
af = @inferred(reinterpret(reshape, Float64, a))
anew = @inferred(reinterpret(reshape, Tuple{Float64,Float64}, vec(af)))
@test anew[1] == a[1]
@test ndims(anew) == 0

# re-reinterpret
a0 = reshape([0x22, 0x44, 0x88, 0xf0, 0x01, 0x02, 0x03, 0x04], 4, 2)
a = reinterpret(reshape, NTuple{4,UInt8}, a0)
@test a == [(0x22, 0x44, 0x88, 0xf0), (0x01, 0x02, 0x03, 0x04)]
@test reinterpret(UInt8, a) == [0x22, 0x44, 0x88, 0xf0, 0x01, 0x02, 0x03, 0x04]
@test reinterpret(reshape, UInt8, a) === a0

# reductions
a = [(1,2,3), (4,5,6)]
ars = reinterpret(reshape, Int, a)
@test sum(ars) == 21
@test sum(ars; dims=1) == [6 15]
@test sum(ars; dims=2) == reshape([5,7,9], (3, 1))
@test sum(ars; dims=(1,2)) == reshape([21], (1, 1))
# also test large sizes for the pairwise algorithm
a = [(k,k+1,k+2) for k = 1:3:4000]
ars = reinterpret(reshape, Int, a)
@test sum(ars) == 8010003

max_stars_repo_name: pengwyn/julia

max_issues_repo_name: pengwyn/julia

max_forks_repo_name: pengwyn/julia


text:
@testset "blockBiCGSTB" begin
	@testset "real matrix" begin
		nrhs = 4;
		n = 100
		A  = sprandn(n,n,.1) + 10*I
		D  = Array(diag(A))
		Af = x ->  A*x
		M  = x ->  D.\x
		rhs = randn(n,nrhs)

		# test flag for early stopping
		xt = blockBiCGSTB(A,rhs,tol=1e-6,maxIter=3,out=2)
		@test xt[2]==-1

		# test handling of zero rhs
		xt = blockBiCGSTB(A,zeros(n,nrhs),tol=1e-6,maxIter=3,out=2)
		@test  all(xt[1].==0)
		@test xt[2]==-9

		x0 = blockBiCGSTB(Matrix(A),rhs,tol=1e-6)
		x1 = blockBiCGSTB(A,rhs,tol=1e-6,out=1)
		x2 = blockBiCGSTB(Af,rhs,tol=1e-6)
		x3 = blockBiCGSTB(Af,rhs,tol=1e-6,maxIter=100,x=randn(size(rhs)))
		x5 = blockBiCGSTB(Af,rhs,tol=1e-6,maxIter=100,M1=M)

		@test norm(A*x0[1]-rhs)/norm(rhs) < 1e-6
		@test norm(A*x1[1]-rhs)/norm(rhs) < 1e-6
		@test norm(A*x2[1]-rhs)/norm(rhs) < 1e-6
		@test norm(A*x3[1]-rhs)/norm(rhs) < 1e-6
		@test norm(A*x5[1]-rhs)/norm(rhs) < 1e-6
		@test norm(x2[1]-x1[1])/norm(x1[1]) < 1e-12
		@test norm(x3[1]-x1[1])/norm(x1[1]) < 1e-5
		@test norm(x0[1]-x1[1])/norm(x1[1]) < 1e-5
	end
	@testset "complex matrix" begin
		nrhs = 4;
		n = 100
		A  = sprandn(n,n,.1) + 10*I + im*(sprandn(n,n,.1) + 10*I )
		D  = Array(diag(A))
		Af = x -> A*x
		M  = x -> D.\x
		rhs = randn(n,nrhs) + 1im * randn(n,nrhs)

		x1 = blockBiCGSTB(A,rhs,tol=1e-6)
		x2 = blockBiCGSTB(Af,rhs,tol=1e-6)
		x3 = blockBiCGSTB(Af,rhs,tol=1e-6,maxIter=200,x=randn(size(rhs))+im*randn(size(rhs)))
		x4 = blockBiCGSTB(Af,rhs,tol=1e-6,maxIter=200,M1=M)

		@test norm(A*x1[1]-rhs)/norm(rhs) < 1e-6
		@test norm(A*x2[1]-rhs)/norm(rhs) < 1e-6
		@test norm(A*x3[1]-rhs)/norm(rhs) < 1e-6
		@test norm(A*x4[1]-rhs)/norm(rhs) < 1e-6
		@test norm(x2[1]-x1[1])/norm(x1[1]) < 1e-12
		@test norm(x3[1]-x1[1])/norm(x1[1]) < 1e-5
	end
end


text:
import Statistics
using Base.Threads
using Memoize

@memoize function get_phenotype_similarity(genotype_1, genotype_2, rel_tolerance, min_samples, max_samples; test_identical = false)
    if genotype_1 == genotype_2 && test_identical == false return 0.0 end

    similarity = estimate(rel_tolerance, min_samples, max_samples) do
	    snapshot = sample_snapshot_id() |> get_snapshot
        sample = get_organisms(snapshot) |> rand
        sample_id = get_id(snapshot, sample)

        sample_similaritiy = (get_fitness(snapshot, sample_id, genotype_1) - get_fitness(snapshot, sample_id, genotype_2))^2

        return sample_similaritiy
    end

    return similarity
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
"""
Created on Wed Apr 15 2020

@author: Yoann Pradat

    CentraleSupelec
    MICS laboratory
    9 rue Juliot Curie, Gif-Sur-Yvette, 91190 France

Implementation of the Maximum-Equalization algorithms and its variants under various constraints and for the β
divergence as presented in C. Fevotte & J. Idier, "Algorithms for nonnegative matrix factorization
with the beta-divergence ", Neural Compuation, 2011.
"""

abstract type AbstractFIParams <: AbstractNMFParams end

"""
    FIParams([,β, scale_W_iter, scale_W_last, α_H, α_W, l₁ratio_H, l₁ratio_W, θ, alg, seed])

Local parameters passed to `nmf_FI`.

# Arguments
- `β::Real`               : parameter of the β-divergence
- `scale_W_iter::Bool`    : scale the components of W to have unit l₁ norm during iterations
- `scale_W_last::Bool`     : scale the components of W to have unit l₁ norm after last iteration
- `α_H::Real`             : constant that multiplies the regularization term on H
- `α_W::Real`             : constant that multiplies the regularization term on W
- `l₁ratio_H::Real`       : l₁/l₂ regularization balance for H (in [0; 1]). 1 is l₁ only, 0 is l₂ only
- `l₁ratio_W::Real`       : l₁/l₂ regularization balance for W (in [0; 1]). 1 is l₁ only, 0 is l₂ only
- `θ::Real`               : θ for mixed update between ME and H. Used only in `nmf_ME`.
- `alg::Symbol`           : symbol to choose the algorithm. Choose
    + :mm for Maximization-Minimization
    + :h  for Heuristic (Maximization-Minimization without exponent)
    + :me for Maximization-Equalization
"""
mutable struct FIParams <: AbstractFIParams
    β::Real
    scale_W_iter::Bool
    scale_W_last::Bool
    α_H::Real
    α_W::Real
    l₁ratio_H::Real
    l₁ratio_W::Real
    θ::Real
    alg::Symbol

    function FIParams(;
        β::Real               = 2,
        scale_W_iter::Bool    = true,
        scale_W_last::Bool    = true,
        α_H::Real             = 0,
        α_W::Real             = 0,
        l₁ratio_H::Real       = 0,
        l₁ratio_W::Real       = 0,
        θ::Real               = 0.95,
        alg::Symbol           = :mm)

        valid_α_H = α_H >= 0
        valid_α_W = α_W >= 0
        valid_l₁ratio_H = l₁ratio_H >= 0 && l₁ratio_H <= 1
        valid_l₁ratio_W = l₁ratio_W >= 0 && l₁ratio_W <= 1
        valid_θ = θ >= 0 && θ <= 1
        valid_alg = alg in [:mm :me :h]

        valid_α_H || throw(ArgumentError("α_H must be > 0"))
        valid_α_W || throw(ArgumentError("α_W must be > 0"))
        valid_l₁ratio_H || throw(ArgumentError("l₁ratio_H must be in [0,1]"))
        valid_l₁ratio_W || throw(ArgumentError("l₁ratio_W must be in [0,1]"))
        valid_θ  || throw(ArgumentError("θ must be between 0 and 1"))
        valid_alg || throw(ArgumentError("alg must be one of :mm, :h or :me"))

        new(β, scale_W_iter, scale_W_last, α_H, α_W, l₁ratio_H, l₁ratio_W, θ, alg)
    end
end

# Overload copy from Base to allow copies FIParams structs
Base.copy(params::FIParams) = FIParams(
        β            = params.β,
        scale_W_iter = params.scale_W_iter,
        scale_W_last = params.scale_W_iter,
        α_W          = params.α_W,
        α_H          = params.α_H,
        l₁ratio_H    = params.l₁ratio_H,
        l₁ratio_W    = params.l₁ratio_W,
        θ            = params.θ,
        alg          = params.alg
    )

"""
    _γ_β(β)

Compute the exponent used in update rules of Maximization-Minimization algorithm.
"""
function _γ_β(β::Real)
    if β < 1
	γ = 1.0 / (2.0-β)
    elseif  β > 2
	γ = 1.0 / (β-1.0)
    else
	γ = 1
    end
    γ
end

"""
    _β_cost(V, W, H, params)

Compute the cost function that is minimized by the NMF alogrithm running with parameters params. It adds up
the β_divergence and regularization terms.

See also: [`β_divergence`](@ref), [`FIParams`](@ref)
"""
function _β_cost(V::Matrix{T}, W::Matrix{T}, H::Matrix{T}, params::FIParams) where T <: Real
    regu_H = params.α_H * (params.l₁ratio_H * norm(H, 1) + (1. - params.l₁ratio_H) * norm(H, 2)) / length(H)
    regu_W = params.α_W * (params.l₁ratio_W * norm(W, 1) + (1. - params.l₁ratio_W) * norm(W, 2)) / length(W)
    β_div  = β_divergence(V, W*H, params.β)

    β_div + regu_H + regu_W
end

"""
    _update_nmf_MM!(V, WH, W, H, params)

Update the H and W matrices in that order using the update rules of Maximization-Minimization algorithm.
"""
function _update_nmf_MM!(V::Matrix{T}, WH::Matrix{T}, W::Matrix{T}, H::Matrix{T}, params::FIParams) where T <: Real
    #### compute regularization
    #### for J_H = ||H||₁, ∇J_H = J_{K,N} with J the matrix filled with ones
    #### for J_H = 0.5*||H||₂, ∇J_H = H
    regu_H = H -> params.α_H * (params.l₁ratio_H * ones(size(H)) + (1. - params.l₁ratio_H) * H) / length(H)
    regu_W = W -> params.α_W * (params.l₁ratio_W * ones(size(W)) + (1. - params.l₁ratio_W) * W) / length(W)

    β = params.β

    #### exponent for multiplicative updates
    γ = _γ_β(β)

    if β > 2
        #### update rules for β > 2
        W .= W .* (max.((V .* WH .^ (β-2)) * transpose(H) - regu_W(W), 0) ./ (WH .^ (β-1) * transpose(H))) .^ γ
        WH .= W * H

        H .= H .* (max.(transpose(W) * (V .* WH .^ (β-2)) - regu_H(H), 0) ./ (transpose(W) * WH .^ (β-1))) .^ γ
        WH .= W * H

    elseif β == 2
        #### update rules for β = 2 can be written more efficiently
        W .= W .* (max.(V * transpose(H) - regu_W(W), 0) ./ (W * (H * transpose(H))))
        H .= H .* (max.((transpose(W) * V - regu_H(H)), 0) ./ ((transpose(W) * W) * H))
        WH .= W * H

    elseif 1 < β < 2
        #### update rules for 1 < β < 2
        #### the update of h_kn is obtained by solving
        #### \sum_f w_{fk} [ (\tilde{v}_f \frac{h_k}{\tilde{h}_k})^{\beta-1} - v_f (\tilde{v}_f \frac{h_k}{\tilde{h}_k})^{\beta-2} ] + \alpha_h \nabla_{h_k} J_k (h_k) = 0 with J_k the component of J_h(h) = \sum_k J_{k}(h_k)

        if params.α_W == 0
            W .= W .* (((V .* WH .^ (β-2)) * transpose(H)) ./ (WH .^ (β-1) * transpose(H)))
        else
            #### solving a X^2 + b X - c = 0 with X = \frac{h_k}{\tilde{h}_k}^{\beta-1}
            a = WH .^ (β-1) * transpose(H)
            b = regu_W(W)
            c = (V .* WH .^ (β-2)) * transpose(H)
            Δ =  b .^ 2 + 4a .* c

            W .= W .* ((- b + sqrt.(Δ)) ./ (2a)) .^ (1. / (β-1))
        end

        WH .= W * H

        if params.α_H == 0
            H .= H .* ((transpose(W) * (V .* WH .^ (β-2))) ./ (transpose(W) * WH .^ (β-1)))
        else
            #### see above
            a = transpose(W) * WH .^ (β-1)
            b = regu_H(H)
            c = transpose(W) * (V .* WH .^ (β-2))
            Δ =  b .^ 2 + 4a .* c

            H .= H .* ((- b + sqrt.(Δ)) ./ (2 .* a)) .^ (1. / (β-1))
        end

        WH .= W * H

    elseif β <= 1
        #### update rules for β <= 1
        W .= W .* (((V .* WH .^ (β-2)) * transpose(H)) ./ (WH .^ (β-1) * transpose(H) + regu_W(W))) .^ γ
        WH .= W * H

        H .= H .* ((transpose(W) * (V .* WH .^ (β-2))) ./ (transpose(W) * WH .^ (β-1) + regu_H(H))) .^ γ
        WH .= W * H
    end

    #### avoid underflow
    W .= max.(W, eps(T))
    H .= max.(H, eps(T))
    WH .= W * H

    if params.scale_W_iter
        scale_col_W!(W, H, 1)
    end
end

"""
    _update_nmf_H!(V, WH, W, H, params)

Update the H and W matrices in that order using the update rules of Heuristic algorithm.
"""
function _update_nmf_H!(V::Matrix{T}, WH::Matrix{T}, W::Matrix{T}, H::Matrix{T}, params::AbstractFIParams) where T <: Real
    #### compute regularization
    #### for J_H = ||H||₁, ∇J_H = J_{K,N} with J the matrix filled with ones
    #### for J_H = 0.5*||H||₂, ∇J_H = H
    regu_H = H -> params.α_H * (params.l₁ratio_H * ones(size(H)) + (1. - params.l₁ratio_H) * H) / length(H)
    regu_W = W -> params.α_W * (params.l₁ratio_W * ones(size(W)) + (1. - params.l₁ratio_W) * W) / length(W)

    if params.β > 2
        #### update rules for β > 2
        W .= W .* (((V .* WH .^ (params.β-2)) * transpose(H) - regu_W(W)) ./ (WH .^ (params.β-1) * transpose(H)))
        WH .= W * H

        H .= H .* ((transpose(W) * (V .* WH .^ (params.β-2)) - regu_H(H)) ./ (transpose(W) * WH .^ (params.β-1)))
        WH .= W * H

    elseif params.β == 2
        #### update rules for β = 2 can be written more efficiently
        W .= W .* ((V * transpose(H) - regu_W(W)) ./ (W * (H * transpose(H))))
        H .= H .* ((transpose(W) * V - regu_H(H)) ./ ((transpose(W) * W) * H))
        WH .= W * H

    elseif 1 < params.β < 2
        #### update rules for 1 < β < 2
        #### the update of h_kn is obtained by solving
        #### \sum_f w_{fk} [ (\tilde{v}_f \frac{h_k}{\tilde{h}_k})^{\beta-1} - v_f (\tilde{v}_f \frac{h_k}{\tilde{h}_k})^{\beta-2} ] + \alpha_h \nabla_{h_k} J_k (h_k) = 0 with J_k the component of J_h(h) = \sum_k J_{k}(h_k)

        if params.α_W == 0
            W .= W .* (((V .* WH .^ (params.β-2)) * transpose(H)) ./ (WH .^ (params.β-1) * transpose(H)))
        else
            #### solving a X^2 + b X - c = 0 with X = \frac{h_k}{\tilde{h}_k}^{\beta-1}
            a = WH .^ (params.β-1) * transpose(H)
            b = regu_W(W)
            c = (V .* WH .^ (params.β-2)) * transpose(H)
            Δ =  b .^ 2 + 4a .* c

            W .= W .* ((- b + sqrt.(Δ)) ./ (2 .* a)) .^ (1. / (params.β-1))
        end

        WH .= W * H

        if params.α_H == 0
            H .= H .* ((transpose(W) * (V .* WH .^ (params.β-2))) ./ (transpose(W) * WH .^ (params.β-1)))
        else
            #### see above
            a = transpose(W) * WH .^ (params.β-1)
            b = regu_H(H)
            c = transpose(W) * (V .* WH .^ (params.β-2))
            Δ =  b .^ 2 + 4a .* c

            H .= H .* ((- b + sqrt.(Δ)) ./ (2a)) .^ (1. / (params.β-1))
        end

        WH .= W * H

    elseif params.β <= 1
        #### update rules for β <= 1
        W .= W .* (((V .* WH .^ (params.β-2)) * transpose(H)) ./ (WH .^ (params.β-1) * transpose(H) + regu_W(W)))
        WH .= W * H

        H .= H .* ((transpose(W) * (V .* WH .^ (params.β-2))) ./ (transpose(W) * WH .^ (params.β-1) + regu_H(H)))
        WH .= W * H
    end

    #### avoid underflow
    W .= max.(W, eps(T))
    H .= max.(H, eps(T))
    WH .= W * H

    if params.scale_W_iter
        scale_col_W!(W, H, 1)
    end
end

"""
    _update_nmf_ME!(V, WH, W, H, params)

Update the H and W matrices in that order using the update rules of Maximization-Equalization algorithm.
"""
function _update_nmf_ME!(V::Matrix{T}, WH::Matrix{T}, W::Matrix{T}, H::Matrix{T}, params::FIParams) where T <: Real
    #### compute regularization
    #### for J_H = ||H||₁, ∇J_H = J_{K,N} with J the matrix filled with ones
    #### for J_H = 0.5*||H||₂, ∇J_H = H
    regu_H = H -> params.α_H * (params.l₁ratio_H * ones(size(H)) + (1. - params.l₁ratio_H) * H) / length(H)
    regu_W = W -> params.α_W * (params.l₁ratio_W * ones(size(W)) + (1. - params.l₁ratio_W) * W) / length(W)

    if params.β == -1
        #### ME always defined
        W_h = W .* (((V .* WH .^ (params.β-2)) * transpose(H)) ./ (WH .^ (params.β-1) * transpose(H) + regu_W(W)))
        W .= ((W_h .^ 2 + 8 * W_h .* W) .^ .5 + W_h) / 4
        WH .= W * H

        H_h = H .* ((transpose(W) * (V .* WH .^ (params.β-2))) ./ (transpose(W) * WH .^ (params.β-1) + regu_H(H)))
        H .= ((H_h.^2 + 8 * H_h .* H) .^ .5 + H_h) / 4
        WH .= W * H

    elseif params.β == 0
        #### ME always defined, equal to H
        W .= W .* (((V .* WH .^ (params.β-2)) * transpose(H)) ./ (WH .^ (params.β-1) * transpose(H) + regu_W(W)))
        WH .= W * H

        H .= H .* ((transpose(W) * (V .* WH .^ (params.β-2))) ./ (transpose(W) * WH .^ (params.β-1) + regu_H(H)))
        WH .= W * H

    elseif params.β == 0.5
        #### ME always defined
        W_h = W .* (((V .* WH .^ (params.β-2)) * transpose(H)) ./ (WH .^ (params.β-1) * transpose(H) + regu_W(W)))
        W .= (sqrt.(W + 8 * W_h) - sqrt.(W)) .^ 2 / 4
        WH .= W * H

        H_h = H .* ((transpose(W) * (V .* WH .^ (params.β-2))) ./ (transpose(W) * WH .^ (params.β-1) + regu_H(H)))
        H .= (sqrt.(H + 8 * H_h) - sqrt.(H)) .^2 / 4
        WH .= W * H

    elseif params.β == 1.5
        #### mixed update prolonged ME and MM

        #### compute W update from Heuristic algorithm
        if params.α_W == 0
            W_h = W .* (((V .* WH .^ (params.β-2)) * transpose(H)) ./ (WH .^ (params.β-1) * transpose(H)))
        else
            a = WH .^ (params.β-1) * transpose(H)
            b = regu_W(W)
            c = (V .* WH .^ (params.β-2)) * transpose(H)
            Δ =  b .^ 2 + 4a .* c

            W_h = W .* ((- b + sqrt.(Δ)) ./ (2a)) .^ (1. / (params.β-1))
        end

        #### compute mixed update
        DW = 12W_h - 3W
        IW = 3W_h .< W
        W .= real.((sqrt.(Complex.(DW)) - sqrt.(W)) .^2 / 4)
        W[IW] .= 0
        W .= W_h + params.θ * (W - W_h)
        WH .= W * H

        #### compute H update from Heuristic algorithm
        if params.α_H == 0
            H_h = H .* ((transpose(W) * (V .* WH .^ (params.β-2))) ./ (transpose(W) * WH .^ (params.β-1)))
        else
            a = transpose(W) * WH .^ (params.β-1)
            b = regu_H(H)
            c = transpose(W) * (V .* WH .^ (params.β-2))
            Δ =  b .^ 2 + 4a .* c

            H_h = H .* ((- b .+ sqrt.(Δ)) ./ (2a)) .^ (1. / (params.β-1))
        end

        #### compute mixed update
        DH = 12H_h - 3H
        IH = 3H_h .< H
        H .= real.((sqrt.(Complex.(DH)) - sqrt.(H)) .^2 / 4)
        H[IH] .= 0
        H .= H_h + params.θ * (H - H_h)
        WH .= W * H

    elseif params.β == 2
        #### mixed update prolonged ME and MM
        W_h = W .* ((V * transpose(H) - regu_W(W)) ./ (W * (H * transpose(H))))
        W .= 2W_h - W
        W[W .<= 0] .= 0
        W .= W_h + params.θ * (W - W_h)
        WH .= W * H

        H_h = H .* ((transpose(W) * V - regu_H(H)) ./ ((transpose(W) * W) * H))
        H .= 2H_h - H
        H[H .<= 0] .= 0
        H .= H_h + params.θ * (H - H_h)
        WH .= W * H

    elseif params.β == 3
        #### mixed update prolonged ME and MM
        W_h = W .* (((V .* WH .^ (params.β-2)) * transpose(H) - regu_W(W)) ./ (WH .^ (params.β-1) * transpose(H)))
        DW = W .* (12W_h - 3W)
        IW = 3W_h .< W
        W .= real.(sqrt.(Complex.(DW)) - W) / 2
        W[IW] .= 0
        W .= W_h + params.θ * (W - W_h)
        WH .= W * H

        H_h = H .* ((transpose(W) * (V .* WH .^ (params.β-2)) - regu_H(H)) ./ (transpose(W) * WH .^ (params.β-1)))

        DH = H .* (12H_h - 3H)
        IH = 3H_h .< H
        H .= real.(sqrt.(Complex.(DH)) - H) / 2
        H[IH] .= 0
        H .= H_h + params.θ * (H - H_h)

        WH .= W * H
    else
        throw(ArgumentError("ME algorithm not implemented for β=$(params.β)"))
    end

    #### avoid underflow
    W .= max.(W, eps(T))
    H .= max.(H, eps(T))
    WH .= W * H

    if params.scale_W_iter
        scale_col_W!(W, H, 1)
    end
end

"""
    _nmf_FI(V, W_init, H_init, global_params, local_params)

Helper function for `nmf_FI`. Runs the FI algorithms with given intial matrices W and H.
"""
function _nmf_FI(V::Matrix{T}, W_init::Matrix{T}, H_init::Matrix{T}, global_params::AbstractNMFParams, local_params::AbstractFIParams) where T <: Real
    W  = W_init[:,:]
    H  = H_init[:,:]
    WH = W * H

    #### metrics and stopping_crit
    func_cost          = (V, W, H) -> _β_cost(V, W, H, local_params)
    func_divg          = (V, W, H) -> β_divergence(V, W * H, local_params.β)
    metrics            = NMFMetrics(V, W_init, H_init, func_cost, func_divg)
    stopping_crit_met  = false
    iteration_is_sane  = true

    #### scale W before first iter if specified
    if local_params.scale_W_iter
        scale_col_W!(W, H, 1)
    end

    while !stopping_crit_met && iteration_is_sane
        #### update WH, W, H
        if local_params.alg == :mm
            _update_nmf_MM!(V, WH, W, H, local_params)
        elseif local_params.alg == :h
            _update_nmf_H!(V, WH, W, H, local_params)
        elseif local_params.alg == :me
            _update_nmf_ME!(V, WH, W, H, local_params)
        end

        #### update metrics
        _update_metrics!(metrics, V, W, H, global_params.stopping_crit)

        #### check sanity of iteration looking for potential divisions by zeros
        iteration_is_sane = _check_sanity_iter(W, H)

        #### check the stopping criterion
        stopping_crit_met = _check_stopping_crit(metrics, global_params)

        #### print intermediate if verbose ≥ 2
        if global_params.verbose ≥ 2 && global_params.max_iter < 10
            @printf("iter %d/%d | cost: %.5g\n", metrics.cur_iter, global_params.max_iter, metrics.cur_cost)
        elseif global_params.verbose ≥ 2 && (metrics.cur_iter % (global_params.max_iter ÷ 10) == 0)
            @printf("iter %d/%d | cost: %.5g\n", metrics.cur_iter, global_params.max_iter, metrics.cur_cost)
        end
    end

    if local_params.scale_W_last
        scale_col_W!(W, H, 1)
    end

    #### convergence information
    if !iteration_is_sane
        @warn "Stopped nmf because sanity check not passed at iteration" metrics.cur_iter
    else
        if global_params.verbose ≥ 1 && metrics.cur_iter == global_params.max_iter && global_params.stopping_crit != :none
            @warn "Maximum number of iterations reached"
        end
        if global_params.verbose ≥ 1 && metrics.cur_iter < global_params.max_iter
            @info "Stopping criterion $(global_params.stopping_crit) met at iter $(metrics.cur_iter)"
        end
    end

    NMFResults{T}(W, H, metrics, global_params, local_params, stopping_crit_met)
end

"""
    nmf_FI(V, global_params, local_params)

Run NMF algorithm using the set of input parameters. For the parameters, see [`Params`](@ref) and [`FIParams`](@ref).

# References
- C. Fevotte & J. Idier, "Algorithms for nonnegative matrix factorization
with the beta-divergence ", Neural Compuation, 2011.
"""
function nmf_FI(V::Matrix{T}, global_params::AbstractNMFParams, local_params::AbstractFIParams) where T <: Real
    # initalization
    W_init, H_init = initialize_nmf(
        V,
        rank = global_params.rank,
        init = global_params.init,
        dist = global_params.dist,
        seed = global_params.seed
    )

    # algorithm
    _nmf_FI(V, W_init, H_init, global_params, local_params)
end


text:
# mf_plots.jl

#==
use this as a template for plotting mfw plots
==#
using Plots, CSV

qValues = [-5:5;]
scales = [3:11;]

filename = "/Users/tehrandavis/Desktop/multifractal_toolbox/exampleTS.txt"

dataframe = CSV.read(filename, header = false)

ts = dataframe.Column1

df = CAPtoolbox.Multifractal.ChhabraJensen(ts, qValues, scales)

p1 = Plots.plot(df.alpha, df.falpha)
p2 = Plots.plot(df.Dq, qValues)

Plots.plot(p1,p2, legend = false)


text:
using Showoff
using Test
using Dates
using Printf

const drops0s = !isdefined(Base, :Ryu)

@testset "Internals" begin
    @test Showoff.@grisu_ccall(1, 2, 3) === nothing
    if isdefined(Showoff, :Grisu)
        @test Showoff.grisu(1.0, Showoff.Grisu.SHORTEST, 2) == (1, 1, false, Showoff.Grisu.DIGITS)
    end

    let x = [1.0, Inf, 2.0, NaN]
        @test Showoff.concrete_minimum(x) == 1.0
        @test Showoff.concrete_maximum(x) == 2.0
    end

    @test_throws ArgumentError Showoff.concrete_minimum([])
    @test_throws ArgumentError Showoff.concrete_maximum([])

    let x = [1.12345, 4.5678]
        @test Showoff.plain_precision_heuristic(x) == 5
        @test Showoff.scientific_precision_heuristic(x) == 6
    end
end

@testset "Formatting" begin
    @test Showoff.format_fixed(-10.0, 0) == "-10"
    @test Showoff.format_fixed(0.012345, 3) == "0.012"
    @test Showoff.format_fixed(Inf, 1) == "∞"
    @test Showoff.format_fixed(-Inf, 1) == "-∞"
    @test Showoff.format_fixed(NaN, 1) == "NaN"
    @test Showoff.format_fixed_scientific(0.0, 1, false) == "0"
    @test Showoff.format_fixed_scientific(Inf, 1, false) == "∞"
    @test Showoff.format_fixed_scientific(-Inf, 1, false) == "-∞"
    @test Showoff.format_fixed_scientific(NaN, 1, false) == "NaN"
    @test Showoff.format_fixed_scientific(0.012345678, 4, true) == "12.346×10⁻³"
    @test Showoff.format_fixed_scientific(0.012345678, 4, false) == "1.2346×10⁻²"
    @test Showoff.format_fixed_scientific(-10.0, 4, false) == (drops0s ? "-1.000×10¹" : "-1.0000×10¹")
    @test Showoff.format_fixed_scientific(-10.0, 4, true) == "-10.000×10⁰"
    @test Showoff.format_fixed_scientific(-10.0, 4, false)[1:end-5] == (drops0s ? @sprintf("%0.4e", -10.0)[1:end-5] : @sprintf("%0.4e", -10.0)[1:end-4])
    @test Showoff.format_fixed_scientific(1.23456e7, 3, false)[1:end-5] == @sprintf("%0.3e", 1.23456e7)[1:end-4]
    @test Showoff.format_fixed_scientific(2.99999999999999956E-16, 2, false) == (drops0s ? "3.0×10⁻¹⁶" : "3.00×10⁻¹⁶")
end

@testset "Showoff" begin
    x = [1.12345, 4.5678]
    @test showoff(x) == ["1.12345", "4.56780"]
    @test showoff([0.0, 50000.0]) == (drops0s ? ["0", "5×10⁴"] : ["0", "5.0×10⁴"])
    @test showoff(x, :plain) == ["1.12345", "4.56780"]
    @test showoff(x, :scientific) == (drops0s ? ["1.12345×10⁰", "4.56780×10⁰"] : ["1.123450×10⁰", "4.567800×10⁰"])
    @test showoff(x, :engineering) == showoff(x, :scientific)
    @test showoff([DateTime("2017-04-11", "yyyy-mm-dd")]) == ["Apr 11, 2017"]
    @test showoff(["a", "b"]) == ["\"a\"", "\"b\""]
    @test showoff([1, 1e39]) == (drops0s ? ["1×10⁰", "1×10³⁹"] : ["1.0×10⁰", "1.0×10³⁹"])
    @test_throws ArgumentError showoff(x, :nevergonnagiveyouup)
    @test showoff([Inf, Inf, NaN]) == ["Inf", "Inf", "NaN"]
end