text:
export NodeSimMeasure, NodeSimScore

"""
Node conservation measure
"""
immutable NodeSimMeasure{T<:AbstractMatrix} <: NetalignMeasure
    S :: T
    function NodeSimMeasure{T}(S::T) where {T<:AbstractMatrix}
        if size(S,1) > size(S,2) error("Matrix dims") end
        new(S)
    end
end
NodeSimMeasure(S::T) where {T<:AbstractMatrix} =  NodeSimMeasure{T}(S)

immutable NodeSimScore <: NetalignScore
    score :: Float64
end

function measure(m::NodeSimMeasure,f::Vector{Int})
    x = 0.0
    S = m.S
    n = size(S,1)
    for i = 1:n
        x += S[i,f[i]]
    end
    NodeSimScore(x/n)
end

dim(m::NodeSimMeasure,d::Int) = size(m.S,d)
dim(m::NodeSimMeasure) = dim(m,2)

NodeSimMeasure(sym::Symbol, args...) =  NodeSimMeasure(Val{sym}(), args...)

"GDV similarity, from the original paper"
function NodeSimMeasure(::Val{:gdvs}, gdv1::AbstractMatrix{T},gdv2::AbstractMatrix{T}) where {T}
    weights = T[1.0, 0.838444533, 0.838444533, 0.838444533, 0.743940642, 0.676889065, 0.743940642, 0.743940642, 0.676889065, 0.743940642, 0.676889065, 0.676889065, 0.676889065, 0.676889065, 0.743940642, 0.676889065, 0.582385175, 0.624879821, 0.676889065, 0.624879821, 0.582385175, 0.582385175, 0.676889065, 0.676889065, 0.676889065, 0.624879821, 0.546456463, 0.676889065, 0.582385175, 0.582385175, 0.546456463, 0.676889065, 0.582385175, 0.582385175, 0.582385175, 0.624879821, 0.582385175, 0.546456463, 0.546456463, 0.624879821, 0.546456463, 0.582385175, 0.546456463, 0.582385175, 0.624879821, 0.624879821, 0.582385175, 0.515333598, 0.546456463, 0.582385175, 0.582385175, 0.515333598, 0.582385175, 0.487881285, 0.624879821, 0.582385175, 0.676889065, 0.582385175, 0.582385175, 0.546456463, 0.515333598, 0.582385175, 0.582385175, 0.515333598, 0.546456463, 0.582385175, 0.546456463, 0.546456463, 0.515333598, 0.624879821, 0.582385175, 0.582385175, 0.676889065]
    if size(gdv1,2) != size(gdv2,2) || size(gdv1,2) != length(weights) error("GDV sizes don't match") end
    @inline kernel(oai,obi,tai,tbi,wi) = wi * abs(oai-obi) / max(tai,tbi)
    D = zeros(T,size(gdv1,1),size(gdv2,1))
    olog1 = log.(gdv1 .+ one(T))
    tlog1 = log.(gdv1 .+ 2one(T));
    olog2 = log.(gdv2 .+ one(T))
    tlog2 = log.(gdv2 .+ 2one(T));
    D = zeros(T,size(gdv1,1),size(gdv2,1))
    for i = 1:length(weights), v = 1:size(gdv2,1), u = 1:size(gdv1,1)
        @inbounds D[u,v] += kernel(olog1[u,i], olog2[v,i], tlog1[u,i], tlog2[v,i], weights[i])
    end
    D ./= sum(weights)
    D .= 1 .- D
    NodeSimMeasure(D)
end

"GDV similarity from the L-GRAAL paper"
function NodeSimMeasure(::Val{:lgraalgdvs}, gdv1::AbstractMatrix{T},gdv2::AbstractMatrix{T}) where {T}
    if size(gdv1,2) != size(gdv2,2) error("GDV sizes don't match") end
    S = zeros(T,size(gdv1,1),size(gdv2,1))
    @inbounds for i = 1:size(gdv1,2), v = 1:size(gdv2,1), u = 1:size(gdv1,1)
        a,b = minmax(gdv1[u,i], gdv2[v,i])
        S[u,v] += ifelse(b>0, a/b, 0.0)
    end
    S ./= size(gdv1,2)
    NodeSimMeasure(S)
end

"The normalized GDV similarity from the H-GRAAL paper"
function NodeSimMeasure(::Val{:hgraalgdvs}, G1::SparseMatrixCSC,G2::SparseMatrixCSC,
                        gdv1::AbstractMatrix{T},gdv2::AbstractMatrix{T},alpha::T=0.5) where {T}
    gm = NodeSimMeasure(:gdvs, gdv1,gdv2)
    deg1 = vec(full(sum(x->T(x),G1,1)))
    deg2 = vec(full(sum(x->T(x),G2,1)))
    maxdeg = maximum(deg1) + maximum(deg2)
    S = gm.S
    @inbounds for v = 1:size(S,2), u = 1:size(S,1)
        S[u,v] = alpha * S[u,v] + (one(T) - alpha) * (deg1[u]+deg2[v]) / maxdeg
    end
    gm
end

type PCA
    mean::AbstractArray
    proj::AbstractMatrix       # projection
    prinvars::AbstractVector # variance of projX*projX' instead of cov(projX) like fit(PCA)
end
"""
Input:
X : m x n matrix
Each column of X is a feature vector
Get the first dimensions containing k fraction of variance

Output:
    PCA object
    with mean, projection matrix, and variances
    U is projection matrix
    projection is T = U' * X
    U*T =~ X, where X is mean centered
"""
function pca(X::AbstractMatrix,k::Real;mu=mean(X,2))
    m,n = size(X)
    if mu!=0
        X = X .- vec(mu)
    else
        mu = zeros(0,0)
    end
    if m <= n # feature size is less than number of samples
        C = X*X'
        res = eigfact!(C)
        d = res.values
        U = res.vectors
        ix = sortperm(d,rev=true)
        d = d[ix]
        clamp!(d, 0.0, Inf)
    else
        C = X'*X
        res = eigfact!(C)
        d = res.values
        V = res.vectors
        ix = sortperm(d,rev=true)
        d = d[ix]
        clamp!(d, 0.0, Inf)
    end
    if !(0.0 <= k <= 1.0) error("k needs to be between 0 and 1") end
    k = findfirst(cumsum(d)./sum(d) .>= k)
    d = d[1:k]
    ix = ix[1:k]
    if m <= n
        U = U[:,ix]
    else
        V = V[:,ix]
        U = X*V # convert evecs from X'*X to X*X'. the evals are the same.
        s = sqrt.(d)
        U = U ./ s'
    end
    PCA(mu, U, d)
end

"The GDV similarity calculated using PCA from Yuriy's dynamic graphlets paper"
function NodeSimMeasure(::Val{:pcagdvs}, gdv1::AbstractMatrix,gdv2::AbstractMatrix)
    X = hcat(gdv1',gdv2') # switch to make columns feature vectors
    X = X[vec(var(X,2)) .>= 1e-5,:] # remove features w/ no variance
    X .-= mean(X,2)
    X ./= std(X,2) # standardize
    res = pca(X,0.99,mu=0) #fit(PCA,X; pratio=0.99, mean=0)
    Y = sqrt.(res.prinvars + 1e-5) .\ res.proj' * X # project and whiten
    R = 1 .- pairwise(CosineDist(), Y[:,1:size(gdv1,1)], Y[:,size(gdv1,1)+1:end]) ./ 2
    NodeSimMeasure(R)
end

"The GDV similarity calculated using CCA"
function NodeSimMeasure(::Val{:ccagdvs}, gdv1::AbstractMatrix,gdv2::AbstractMatrix)
    X = gdv1'
    Y = gdv2'
    res = fit(CCA,X,Y)
    R = 1 .- pairwise(CosineDist(), xtransform(res,X), ytransform(res,Y)) ./ 2
    NodeSimMeasure(R)
end

"Relative similarity between v1 and v2: 1 - abs(a-b)/max(a,b)"
function NodeSimMeasure(::Val{:relative}, v1::AbstractVector,v2::AbstractVector)
    m = length(v1)
    n = length(v2)
    S = zeros(Float64,m,n)
    for i = 1:m, j = 1:n
        S[i,j] = 1 - abs(v1[i] - v2[j])/max(v1[i],v2[j])
    end
    NodeSimMeasure(S)
end

"Degree similarity"
NodeSimMeasure(::Val{:degree}, G1::SparseMatrixCSC,G2::SparseMatrixCSC) =
    NodeSimMeasure(:relative, vec(sum(G1,1)), vec(sum(G2,1)))

"Relative clustering coefficient"
NodeSimMeasure(::Val{:clustering}, G1::SparseMatrixCSC,G2::SparseMatrixCSC) =
    NodeSimMeasure(:relative, clustercoeffs(G1), clustercoeffs(G2))

"Eccentricity similarity"
function NodeSimMeasure(::Val{:eccentricity}, G1::SparseMatrixCSC,
                        G2::SparseMatrixCSC)
    D1 = hcat(map(i->dijkstra(G1,i)[1],1:size(G1,1))...)
    D2 = hcat(map(i->dijkstra(G2,i)[1],1:size(G2,1))...)
    NodeSimMeasure(:relative, maximum(D1,1), maximum(D2,1))
end

"Performs the migraal confidence transform"
function NodeSimMeasure(::Val{:migraal}, Xs::AbstractMatrix...)
    C = zeros(Float64,size(Xs[1]))
    for X in Xs
        for i = 1:size(X,1)
            for j = 1:size(X,2)
                sless = 0
                for k = 1:size(X,2)
                    sless += Int(X[i,k] < X[i,j])
                end
                C[i,j] = sless/size(X,2)
            end
        end
    end
    C ./= length(Xs)
    NodeSimMeasure(C)
end

"""Transforms blast e-values to nodes similarities.
# Arguments
- `clamp = 1e-101` : Clamp E-values to a minimum of `clamp`."""
function NodeSimMeasure(::Val{:evalues}, E::AbstractMatrix, clamp=1e-101)
    # -log(e-value) transform
    B = -log.(E)
    # Clamp E-value to a minimum of `clamp`
    # That is, clamp -log(E-value) to max of -log(clampval)
    # Otherwise we get issues with -log(E-value) being Inf or some
    # ridiculously huge number
    clamp!(B, 0.0, -log(clamp))
    # Normalize by dividing by maximum value
    B ./= maximum(B)
    NodeSimMeasure(B)
end

"Transforms blast e-values to node similarities"
function NodeSimMeasure(::Val{:evalues}, E::SparseMatrixCSC, clamp=1e-101)
    B = copy(E)
    # Assume structural sparsity, ie if B[i,j] is not in the structural non-zeros
    # of B, then there is no given E-value between i and j
    # But if B[i,j] is, then take the -log(E-value) and clamp the result
    @. B.nzval = -log(B.nzval)
    clamp!(B.nzval, 0.0, -log(clamp))
    B.nzval ./= maximum(B.nzval)
    NodeSimMeasure(B)
end

"GHOST node similiarity"
function NodeSimMeasure(::Val{:ghost}, G1::SparseMatrixCSC,
                        G2::SparseMatrixCSC,k::Integer=4;verbose=true)
    m = size(G1,1)
    n = size(G2,1)
    gmeas = GhostMeasure(G1,G2,k,verbose=verbose)
    verbose && println("Calculating GHOST node similarities")
    R = zeros(Float64,m,n)
    for v = 1:n, u = 1:m
        R[u,v] = score(meas,u,v)
    end
    NodeSimMeasure(R)
end

# interaction matrix = min(G1.expci[i],G2.expci[j])/max(G1.maxdeg,G2.maxdeg)
function NodeSimMeasure(::Val{:interactionscore}, G1::SparseMatrixCSC,G2::SparseMatrixCSC)
    NodeSimMeasure(expci(G1,G2))
end


text:
using StaticArrays, HybridArrays, Test, LinearAlgebra

@testset "AbstractArray interface" begin
    @testset "size and length" begin
        M = HybridMatrix{2, StaticArrays.Dynamic()}([1 2 3; 4 5 6])

        @test length(M) == 6
        @test size(M) == (2, 3)
        @test Base.isassigned(M, 2, 2) == true
    end

    @testset "reshape" begin
        # TODO?
    end

    @testset "convertsions" begin
        M = HybridMatrix{2, StaticArrays.Dynamic()}([1 2 3; 4 5 6])
        @test convert(Array, M) === M.data
        @test_broken convert(Matrix, M) === M.data
        @test convert(Matrix{Int}, M) === M.data

        @test Array(M) == M
        @test Array(M) !== M.data
        @test Matrix(M) == M
        @test Matrix(M) !== M.data
        @test Matrix{Int}(M) == M
        @test Matrix{Int}(M) !== M.data
        @test_throws MethodError Vector(M)
    end

    @testset "copy" begin
        M = HybridMatrix{2, StaticArrays.Dynamic()}([1 2; 3 4])
        @test @inferred(copy(M))::HybridMatrix == M
        @test copy(M).data !== M.data
    end

    @testset "similar" begin
        M = HybridMatrix{2, StaticArrays.Dynamic(), Int}([1 2; 3 4])

        @test isa(@inferred(similar(M)), HybridMatrix{2, StaticArrays.Dynamic(), Int})
        @test isa(@inferred(similar(M, Float64)), HybridMatrix{2, StaticArrays.Dynamic(), Float64})
    end

end

ext: jl


text:
using DrWatson

DrWatson.greet()


text:
@testset "HLADataset(Rousseau)" begin
    @test HLADataset("Rousseau") isa HLADataset

    dataset = HLADataset("Rousseau")
    @test dataset.name == "Rousseau"
    @test dataset.data[1] isa HLAData
    @test dataset.data[2] isa HLAData
    @test dataset isa AbstractHLADataset
end


text:
using .Distances

@adjoint function sqeuclidean(x::AbstractVector, y::AbstractVector)
  δ = x .- y
  return sum(abs2, δ), function(Δ::Real)
    x̄ = (2 * Δ) .* δ
    return x̄, -x̄
  end
end

@adjoint function colwise(s::SqEuclidean, x::AbstractMatrix, y::AbstractMatrix)
  return colwise(s, x, y), function (Δ::AbstractVector)
    x̄ = 2 .* Δ' .* (x .- y)
    return nothing, x̄, -x̄
  end
end

@adjoint function pairwise(s::SqEuclidean, x::AbstractMatrix, y::AbstractMatrix)
  return pairwise(s, x, y), function(Δ)
    x̄ = 2 .* (x * Diagonal(vec(sum(Δ; dims=2))) .- y * Δ')
    ȳ = 2 .* (y * Diagonal(vec(sum(Δ; dims=1))) .- x * Δ)
    return nothing, x̄, ȳ
  end
end

@adjoint function pairwise(s::SqEuclidean, X::AbstractMatrix)
  D = pairwise(s, X)
  return D, function(Δ)
    d1, d2 = Diagonal(vec(sum(Δ; dims=1))), Diagonal(vec(sum(Δ; dims=2)))
    return (nothing, X * (2 .* (d1 .+ d2 .- Δ .- Δ')))
  end
end


text:
using CLIMA
using CLIMA.Mesh.Topologies: StackedCubedSphereTopology, cubedshellwarp, grid1d
using CLIMA.Mesh.Grids
using CLIMA.Mesh.Filters
using CLIMA.DGmethods: DGModel, init_ode_state
using CLIMA.DGmethods.NumericalFluxes: Rusanov, CentralGradPenalty,
                                       CentralNumericalFluxDiffusive
using CLIMA.ODESolvers: solve!, gettime
using CLIMA.LowStorageRungeKuttaMethod: LSRK144NiegemannDiehlBusch
using CLIMA.VTK: writevtk, writepvtu
using CLIMA.GenericCallbacks: EveryXWallTimeSeconds, EveryXSimulationSteps
using CLIMA.MPIStateArrays: euclidean_distance
using CLIMA.PlanetParameters: R_d, grav, MSLP, planet_radius, cp_d, cv_d, day
using CLIMA.MoistThermodynamics: air_density, total_energy, soundspeed_air, internal_energy, air_temperature
using CLIMA.Atmos: AtmosModel, SphericalOrientation, NoReferenceState,
                   DryModel, NoRadiation, NoSubsidence, NoFluxBC,
                   ConstantViscosityWithDivergence,
                   vars_state, vars_aux,
                   Gravity, Coriolis,
                   HydrostaticState, IsothermalProfile
using CLIMA.VariableTemplates: flattenednames

using MPI, Logging, StaticArrays, LinearAlgebra, Printf, Dates, Test

const output_vtk = true
const ArrayType = CLIMA.array_type()

function main()
  CLIMA.init()
  mpicomm = MPI.COMM_WORLD

  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = Dict("DEBUG" => Logging.Debug,
                  "WARN"  => Logging.Warn,
                  "ERROR" => Logging.Error,
                  "INFO"  => Logging.Info)[ll]

  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))

  polynomialorder = 5
  numelem_horz = 6
  numelem_vert = 8
  timeend = 60 # 400day
  outputtime = 2day

  for FT in (Float64,)

    run(mpicomm, polynomialorder, numelem_horz, numelem_vert,
        timeend, outputtime, ArrayType, FT)
  end
end

function run(mpicomm, polynomialorder, numelem_horz, numelem_vert,
             timeend, outputtime, ArrayType, FT)

  setup = HeldSuarezSetup{FT}()

  vert_range = grid1d(FT(planet_radius), FT(planet_radius + setup.domain_height), nelem = numelem_vert)
  topology = StackedCubedSphereTopology(mpicomm, numelem_horz, vert_range)

  grid = DiscontinuousSpectralElementGrid(topology,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = polynomialorder,
                                          meshwarp = cubedshellwarp)

  model = AtmosModel(SphericalOrientation(),
                     HydrostaticState(IsothermalProfile(setup.T_initial), FT(0)),
                     ConstantViscosityWithDivergence(FT(0)),
                     DryModel(),
                     NoRadiation(),
                     NoSubsidence{FT}(),
                     (Gravity(), Coriolis(), held_suarez_forcing!),
                     NoFluxBC(),
                     setup)

  dg = DGModel(model, grid, Rusanov(),
               CentralNumericalFluxDiffusive(), CentralGradPenalty())

  # determine the time step
  element_size = (setup.domain_height / numelem_vert)
  acoustic_speed = soundspeed_air(FT(315))
  lucas_magic_factor = 5
  dt = lucas_magic_factor * min_node_distance(grid) / acoustic_speed

  Q = init_ode_state(dg, FT(0))
  lsrk = LSRK144NiegemannDiehlBusch(dg, Q; dt = dt, t0 = 0)

  filterorder = 14
  filter = ExponentialFilter(grid, 0, filterorder)
  cbfilter = EveryXSimulationSteps(1) do
    Filters.apply!(Q, 1:size(Q, 2), grid, filter)
    nothing
  end

  eng0 = norm(Q)
  @info @sprintf """Starting
                    ArrayType       = %s
                    FT              = %s
                    polynomialorder = %d
                    numelem_horz    = %d
                    numelem_vert    = %d
                    filterorder     = %d
                    dt              = %.16e
                    norm(Q₀)        = %.16e
                    """ "$ArrayType" "$FT" polynomialorder numelem_horz numelem_vert filterorder dt eng0

  # Set up the information callback
  starttime = Ref(now())
  cbinfo = EveryXWallTimeSeconds(60, mpicomm) do (s=false)
    if s
      starttime[] = now()
    else
      energy = norm(Q)
      runtime = Dates.format(convert(DateTime, now() - starttime[]), dateformat"HH:MM:SS")
      @info @sprintf """Update
                        simtime = %.16e
                        runtime = %s
                        norm(Q) = %.16e
                        """ gettime(lsrk) runtime energy
    end
  end
  callbacks = (cbinfo, cbfilter)

  if output_vtk
    # create vtk dir
    vtkdir = "vtk_heldsuarez" *
      "_poly$(polynomialorder)_horz$(numelem_horz)_vert$(numelem_vert)" *
      "_filter$(filterorder)_$(ArrayType)_$(FT)"
    mkpath(vtkdir)

    vtkstep = 0
    # output initial step
    do_output(mpicomm, vtkdir, vtkstep, dg, Q, model)

    # setup the output callback
    cbvtk = EveryXSimulationSteps(floor(outputtime / dt)) do
      vtkstep += 1
      Qe = init_ode_state(dg, gettime(lsrk))
      do_output(mpicomm, vtkdir, vtkstep, dg, Q, model)
    end
    callbacks = (callbacks..., cbvtk)
  end

  solve!(Q, lsrk; timeend=timeend, callbacks=callbacks)

  # final statistics
  engf = norm(Q)
  @info @sprintf """Finished
  norm(Q)                 = %.16e
  norm(Q) / norm(Q₀)      = %.16e
  norm(Q) - norm(Q₀)      = %.16e
  """ engf engf/eng0 engf-eng0
end

Base.@kwdef struct HeldSuarezSetup{FT}
  p_ground::FT = MSLP
  T_initial::FT = 255
  domain_height::FT = 30e3
end

function (setup::HeldSuarezSetup)(state, aux, coords, t)
  # callable to set initial conditions
  FT = eltype(state)

  r = norm(coords, 2)
  h = r - FT(planet_radius)

  scale_height = R_d * setup.T_initial / grav
  p = setup.p_ground * exp(-h / scale_height)

  state.ρ = air_density(setup.T_initial, p)
  state.ρu = SVector{3, FT}(0, 0, 0)
  state.ρe = state.ρ * (internal_energy(setup.T_initial) + aux.orientation.Φ)
  nothing
end

function held_suarez_forcing!(source, state, aux, t::Real)
  FT = eltype(state)

  ρ = state.ρ
  ρu = state.ρu
  ρe = state.ρe
  coord = aux.coord
  Φ = aux.orientation.Φ
  e = ρe / ρ
  u = ρu / ρ
  e_int = e - u' * u / 2 - Φ
  T = air_temperature(e_int)
  # Held-Suarez constants
  k_a = FT(1 / (40 * day))
  k_f = FT(1 / day)
  k_s = FT(1 / (4 * day))
  ΔT_y = FT(60)
  Δθ_z = FT(10)
  T_equator = FT(315)
  T_min = FT(200)
  σ_b = FT(7 / 10)
  r = norm(coord, 2)
  @inbounds λ = atan(coord[2], coord[1])
  @inbounds φ = asin(coord[3] / r)
  h = r - FT(planet_radius)
  scale_height = FT(7000) #from Smolarkiewicz JAS 2001 paper
  σ = exp(-h / scale_height)
  # TODO: use
  #  p = air_pressure(T, ρ)
  #  σ = p/p0
  exner_p = σ ^ (R_d / cp_d)
  Δσ = (σ - σ_b) / (1 - σ_b)
  height_factor = max(0, Δσ)
  T_equil = (T_equator - ΔT_y * sin(φ) ^ 2 - Δθ_z * log(σ) * cos(φ) ^ 2 ) * exner_p
  T_equil = max(T_min, T_equil)
  k_T = k_a + (k_s - k_a) * height_factor * cos(φ) ^ 4
  k_v = k_f * height_factor
  source.ρu += -k_v * ρu
  source.ρe += -k_T * ρ * cv_d * (T - T_equil) - k_v * ρu' * ρu / ρ
end

function do_output(mpicomm, vtkdir, vtkstep, dg, Q, model, testname = "heldsuarez")
  ## name of the file that this MPI rank will write
  filename = @sprintf("%s/%s_mpirank%04d_step%04d",
                      vtkdir, testname, MPI.Comm_rank(mpicomm), vtkstep)

  statenames = flattenednames(vars_state(model, eltype(Q)))
  auxnames = flattenednames(vars_aux(model, eltype(Q)))
  writevtk(filename, Q, dg, statenames, dg.auxstate, auxnames)

  ## Generate the pvtu file for these vtk files
  if MPI.Comm_rank(mpicomm) == 0
    ## name of the pvtu file
    pvtuprefix = @sprintf("%s/%s_step%04d", vtkdir, testname, vtkstep)

    ## name of each of the ranks vtk files
    prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
      @sprintf("%s_mpirank%04d_step%04d", testname, i - 1, vtkstep)
    end

    writepvtu(pvtuprefix, prefixes, (statenames..., auxnames...))

    @info "Done writing VTK: $pvtuprefix"
  end
end

main()

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
## ------------------------------------------------------------------
# References

## ------------------------------------------------------------------
# Meta

refdict(r::RBRef) = getproperty(r, :dict, nothing)
refdict!(r::RBRef, dict::AbstractDict) = setproperty!(r, :dict, dict)

get_bibkey(r::RBRef) = getproperty(r, :bibkey, "")
set_bibkey!(r::RBRef, key::String) = setproperty!(r, :bibkey, strip(key))

get_author(r::RBRef) = getproperty(r, :author, "")
set_author!(r::RBRef, author::String) = setproperty!(r, :author, strip(author))

get_year(r::RBRef) = getproperty(r, :year, "")
set_year!(r::RBRef, year::String) = setproperty!(r, :year, strip(year))

get_title(r::RBRef) = getproperty(r, :title, "")
set_title!(r::RBRef, title::String) = setproperty!(r, :title, strip(title))

get_doi(r::RBRef) = getproperty(r, :doi, "")
set_doi!(r::RBRef, doi::String) = setproperty!(r, :doi, _doi_to_url(doi))

## ------------------------------------------------------------------
# Base
function Base.show(io::IO, r::RBRef)
    print(io, "RBRef")
    for meta in [:title, :bibkey, :author, :year, :doi]
        val = getproperty(r, meta)
        isempty(val) && continue
        print(io, "\n - ", meta, ": \"", _preview(io, val), "\"")
    end
    print(io, "\n", _preview(io, "-"^70))
end


text:
using Distributions, CSV, DataFrames, StatsModels
using Random, LinearAlgebra, SparseArrays

include("testutils.jl")

# generate test data
n = 200
p = 5
d = 4
ncats = 3

Random.seed!(33)
ctotal = 1 .+ rand(Poisson(5),n)
vs = rand(n,p-ncats+1)
covarsdf = DataFrame(vs, [:y, :x, :z])
covarsdf[:cat] = CategoricalArray(rand(["$i" for i=1:ncats], n))
ηfn(vi,g) = exp.([0 + i*sum(vi) - parse(Float64, g) for i=1:d])
q = [ηfn(vs[i,:], covarsdf[i,:cat]) for i=1:n]
for i=1:n
  q[i] ./= sum(q[i])
end

counts = convert(SparseMatrixCSC{Float64,Int}, hcat(broadcast((qi,mi)->rand(Multinomial(mi, qi)),q,ctotal)...)')
countsint = convert(Matrix{Int},counts) # used for testing the int eltype case

CSV.write(joinpath(testdir,"data","countsint.csv.gz"),DataFrame(countsint))
CSV.write(joinpath(testdir,"data","covarsdf.csv.gz"),covarsdf)

# construct equivalent covars matrix so we can show how that api works too
covars = ModelMatrix(ModelFrame(@formula(y ~ x + z + cat + y), covarsdf)).m[:,2:end]

# used for testing predict
newcovars = covars[1:10,:]

# used for testing srproj and fit(CIR)
projdir = size(covars,2)

γdistrom = 1.0

# generate R benchmark
using RCall
R"library(textir)"
R"library(Matrix)"
R"cl <- makeCluster(2,type=\"FORK\")"
R"fits <- dmr(cl, $covars, $counts, gamma=$γdistrom, verb=0)"
R"stopCluster(cl)"
coefsRdistrom = rcopy(R"as.matrix(coef(fits))")
zRdistrom = rcopy(R"as.matrix(srproj(fits,$counts))")
z1Rdistrom = rcopy(R"as.matrix(srproj(fits,$counts,$projdir))")
predictRdistrom = rcopy(R"as.matrix(predict(fits,$newcovars,type=\"response\"))")

# save R benchmark
CSV.write(joinpath(testdir,"data","dmr_coefsRdistrom.csv.gz"),DataFrame(Matrix(coefsRdistrom)))
CSV.write(joinpath(testdir,"data","dmr_zRdistrom.csv.gz"),DataFrame(zRdistrom))
CSV.write(joinpath(testdir,"data","dmr_z1Rdistrom.csv.gz"),DataFrame(z1Rdistrom))
CSV.write(joinpath(testdir,"data","dmr_predictRdistrom.csv.gz"),DataFrame(predictRdistrom))


text:
@testset "Writing" begin
    @testset "Variables" begin
        m = CP.FlatZinc.Model()
        @test MOI.is_empty(m)

        b, b_c = MOI.add_constrained_variable(m, MOI.LessThan(0.0))
        c, c_c = MOI.add_constrained_variable(m, MOI.EqualTo(0.0))
        d, d_c = MOI.add_constrained_variable(m, MOI.Interval(0.0, 1.0))
        f, f_c = MOI.add_constrained_variable(m, MOI.LessThan(0))
        g, g_c = MOI.add_constrained_variable(m, MOI.EqualTo(0))
        h, h_c = MOI.add_constrained_variable(m, MOI.Interval(0, 1))
        i, i_c = MOI.add_constrained_variable(m, MOI.EqualTo(false))
        j = MOI.add_variable(m)

        @test MOI.is_valid(m, b)
        @test MOI.is_valid(m, b_c)
        @test MOI.is_valid(m, c)
        @test MOI.is_valid(m, c_c)
        @test MOI.is_valid(m, d)
        @test MOI.is_valid(m, d_c)
        @test MOI.is_valid(m, f)
        @test MOI.is_valid(m, f_c)
        @test MOI.is_valid(m, g)
        @test MOI.is_valid(m, g_c)
        @test MOI.is_valid(m, h)
        @test MOI.is_valid(m, h_c)
        @test MOI.is_valid(m, i)
        @test MOI.is_valid(m, i_c)
        @test MOI.is_valid(m, j)

        # Generate the FZN file.
        io = IOBuffer(truncate=true)
        write(io, m)
        fzn = String(take!(io))

        @test fzn == """var -1.7976931348623157e308..0.0: x1 :: output_var;
            var float: x2 :: output_var = 0.0;
            var 0.0..1.0: x3 :: output_var;
            var -9223372036854775808..0: x4 :: output_var;
            var int: x5 :: output_var = 0;
            var 0..1: x6 :: output_var;
            var bool: x7 :: output_var = false;
            var float: x8 :: output_var;

            solve satisfy;
            """
    end

    @testset "Constraints: CP.MinimumAmong / CP.MaximumAmong" begin
        m = CP.FlatZinc.Model()
        @test MOI.is_empty(m)

        # Create variables.
        x, x_int =
            MOI.add_constrained_variables(m, [MOI.Integer() for _ in 1:5])
        y = MOI.add_variables(m, 5)

        @test !MOI.is_empty(m)
        for i in 1:5
            @test MOI.is_valid(m, x[i])
            @test MOI.is_valid(m, x_int[i])
            @test MOI.is_valid(m, y[i])
        end

        # Don't set names to check whether they are made unique before
        # generating the model.

        # Add constraints.
        c1 = MOI.add_constraint(
            m,
            MOI.VectorOfVariables([x[1], x[2], x[3]]),
            CP.MaximumAmong(2),
        )
        c2 = MOI.add_constraint(
            m,
            MOI.VectorOfVariables(x),
            CP.MinimumAmong(4),
        )
        c3 = MOI.add_constraint(
            m,
            MOI.VectorOfVariables([y[1], y[2], y[3]]),
            CP.MaximumAmong(2),
        )
        c4 = MOI.add_constraint(
            m,
            MOI.VectorOfVariables(y),
            CP.MinimumAmong(4),
        )

        @test MOI.is_valid(m, c1)
        @test MOI.is_valid(m, c2)
        @test MOI.is_valid(m, c3)
        @test MOI.is_valid(m, c4)

        # Test some attributes for these constraints.
        @test MOI.get(m, MOI.ConstraintFunction(), c1) ==
              MOI.VectorOfVariables([x[1], x[2], x[3]])
        @test MOI.get(m, MOI.ConstraintSet(), c1) == CP.MaximumAmong(2)
        @test MOI.get(
            m,
            MOI.ListOfConstraintIndices{
                MOI.VectorOfVariables,
                CP.Element{Int},
            }(),
        ) == MOI.ConstraintIndex[]
        @test MOI.get(
            m,
            MOI.ListOfConstraintIndices{
                MOI.VectorOfVariables,
                CP.MaximumAmong,
            }(),
        ) == [c1, c3]
        @test length(MOI.get(m, MOI.ListOfConstraintTypesPresent())) == 3

        # Generate the FZN file.
        io = IOBuffer(truncate=true)
        write(io, m)
        fzn = String(take!(io))

        @test fzn == """var int: x1 :: output_var;
            var int: x2 :: output_var;
            var int: x3 :: output_var;
            var int: x4 :: output_var;
            var int: x5 :: output_var;
            var float: x6 :: output_var;
            var float: x7 :: output_var;
            var float: x8 :: output_var;
            var float: x9 :: output_var;
            var float: x10 :: output_var;

            constraint array_int_maximum(x1, [x2, x3]);
            constraint array_int_minimum(x1, [x2, x3, x4, x5]);
            constraint array_float_maximum(x6, [x7, x8]);
            constraint array_float_minimum(x6, [x7, x8, x9, x10]);

            solve satisfy;
            """

        # Test that the names have been correctly transformed.
        for v in x
            vn = MOI.get(m, MOI.VariableName(), v)
            @test match(r"^x\d+$", vn) !== nothing
        end
    end

    @testset "Constraints: CP.Element" begin
        m = CP.FlatZinc.Model()
        @test MOI.is_empty(m)

        # Create variables.
        x, x_int =
            MOI.add_constrained_variables(m, [MOI.Integer() for _ in 1:2])
        y, y_bool =
            MOI.add_constrained_variables(m, [MOI.ZeroOne() for _ in 1:2])
        z = MOI.add_variables(m, 2)

        @test !MOI.is_empty(m)
        for i in 1:2
            @test MOI.is_valid(m, x[i])
            @test MOI.is_valid(m, x_int[i])
            @test MOI.is_valid(m, y[i])
            @test MOI.is_valid(m, y_bool[i])
            @test MOI.is_valid(m, z[i])
        end

        # Don't set names to check whether they are made unique before
        # generating the model.

        # Add constraints.
        c1 = MOI.add_constraint(
            m,
            MOI.VectorOfVariables(x),
            CP.Element([6, 5, 4]),
        )
        c2 = MOI.add_constraint(
            m,
            MOI.VectorOfVariables(y),
            CP.Element([true, false]),
        )
        c3 = MOI.add_constraint(
            m,
            MOI.VectorOfVariables(z),
            CP.Element([1.0, 2.0]),
        )

        @test MOI.is_valid(m, c1)
        @test MOI.is_valid(m, c2)
        @test MOI.is_valid(m, c3)

        # Test some attributes for these constraints.
        @test MOI.get(m, MOI.ConstraintFunction(), c1) ==
              MOI.VectorOfVariables(x)
        @test MOI.get(m, MOI.ConstraintSet(), c1) == CP.Element([6, 5, 4])
        @test MOI.get(
            m,
            MOI.ListOfConstraintIndices{
                MOI.VectorOfVariables,
                CP.Element{Int},
            }(),
        ) == [c1]
        @test length(MOI.get(m, MOI.ListOfConstraintTypesPresent())) == 5

        # Generate the FZN file.
        io = IOBuffer(truncate=true)
        write(io, m)
        fzn = String(take!(io))

        @test fzn == """var int: x1 :: output_var;
            var int: x2 :: output_var;
            var bool: x3 :: output_var;
            var bool: x4 :: output_var;
            var float: x5 :: output_var;
            var float: x6 :: output_var;

            array [1..3] of int: ARRAY0 = [6, 5, 4];
            array [1..2] of bool: ARRAY1 = [true, false];
            array [1..2] of float: ARRAY2 = [1.0, 2.0];

            constraint array_int_element(x2, ARRAY0, x1);
            constraint array_bool_element(x4, ARRAY1, x3);
            constraint array_float_element(x6, ARRAY2, x5);

            solve satisfy;
            """

        # Test that the names have been correctly transformed.
        for v in [x..., y..., z...]
            vn = MOI.get(m, MOI.VariableName(), v)
            @test match(r"^x\d+$", vn) !== nothing
        end
    end

    @testset "Constraints: CP.DifferentFrom" begin
        m = CP.FlatZinc.Model()
        @test MOI.is_empty(m)

        # Create variables.
        x, x_int =
            MOI.add_constrained_variables(m, [MOI.Integer() for _ in 1:2])
        y = MOI.add_variables(m, 2)

        @test !MOI.is_empty(m)
        for i in 1:2
            @test MOI.is_valid(m, x[i])
            @test MOI.is_valid(m, x_int[i])
            @test MOI.is_valid(m, y[i])
        end

        # Don't set names to check whether they are made unique before
        # generating the model.

        # Add constraints.
        c1 = MOI.add_constraint(
            m,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, 1], x), 0),
            CP.DifferentFrom(2),
        )
        c2 = MOI.add_constraint(
            m,
            MOI.ScalarAffineFunction(
                MOI.ScalarAffineTerm.([1.0, 2.0], y),
                0.0,
            ),
            CP.DifferentFrom(2.0),
        )

        @test MOI.is_valid(m, c1)
        @test MOI.is_valid(m, c2)

        # Test some attributes for these constraints.
        @test length(MOI.get(m, MOI.ListOfConstraintTypesPresent())) == 3

        # Generate the FZN file.
        io = IOBuffer(truncate=true)
        write(io, m)
        fzn = String(take!(io))

        @test fzn == """var int: x1 :: output_var;
            var int: x2 :: output_var;
            var float: x3 :: output_var;
            var float: x4 :: output_var;

            constraint int_lin_ne([1, 1], [x1, x2], 2);
            constraint float_lin_ne([1.0, 2.0], [x3, x4], 2.0);

            solve satisfy;
            """

        # Test that the names have been correctly transformed.
        for v in [x..., y...]
            vn = MOI.get(m, MOI.VariableName(), v)
            @test match(r"^x\d+$", vn) !== nothing
        end
    end

    @testset "Constraints: CP.Domain" begin
        m = CP.FlatZinc.Model()
        @test MOI.is_empty(m)

        # Create variable.
        x, x_int = MOI.add_constrained_variable(m, MOI.Integer())
        @test !MOI.is_empty(m)
        @test MOI.is_valid(m, x)
        @test MOI.is_valid(m, x_int)

        # Don't set names to check whether they are made unique before
        # generating the model.

        # Create constraint.
        c = MOI.add_constraint(m, x, CP.Domain(Set([0, 1, 2])))
        @test MOI.is_valid(m, c)

        # Test some attributes for this constraint.
        @test length(MOI.get(m, MOI.ListOfConstraintTypesPresent())) == 2

        # Generate the FZN file.
        io = IOBuffer(truncate=true)
        write(io, m)
        fzn = String(take!(io))

        @test fzn == """var int: x1 :: output_var;

            set of int: SET0 = {0, 2, 1};

            constraint set_in(x1, SET0);

            solve satisfy;
            """

        # Test that the name has been correctly transformed.
        xn = MOI.get(m, MOI.VariableName(), x)
        @test match(r"^x\d+$", xn) !== nothing
    end

    @testset "Constraints: CP.Interval" begin
        m = CP.FlatZinc.Model()
        @test MOI.is_empty(m)

        # Create variable.
        x = MOI.add_variable(m)
        @test !MOI.is_empty(m)
        @test MOI.is_valid(m, x)

        # Don't set names to check whether they are made unique before
        # generating the model.

        # Create constraint.
        c = MOI.add_constraint(m, x, MOI.Interval(1.0, 2.0))
        @test MOI.is_valid(m, c)

        # Test some attributes for this constraint.
        @test length(MOI.get(m, MOI.ListOfConstraintTypesPresent())) == 1

        # Generate the FZN file.
        io = IOBuffer(truncate=true)
        write(io, m)
        fzn = String(take!(io))

        @test fzn == """var float: x1 :: output_var;

            constraint float_in(x1, 1.0, 2.0);

            solve satisfy;
            """

        # Test that the name has been correctly transformed.
        xn = MOI.get(m, MOI.VariableName(), x)
        @test match(r"^x\d+$", xn) !== nothing
    end

    @testset "Constraints: MOI.ScalarAffineFunction in MOI.EqualTo / MOI.LessThan / CP.Strictly{MOI.LessThan} / CP.DifferentFrom" begin
        m = CP.FlatZinc.Model()
        @test MOI.is_empty(m)

        # Create variables.
        x, x_int =
            MOI.add_constrained_variables(m, [MOI.Integer() for _ in 1:2])
        y, y_bool =
            MOI.add_constrained_variables(m, [MOI.ZeroOne() for _ in 1:2])
        z = MOI.add_variables(m, 2)

        @test !MOI.is_empty(m)
        for i in 1:2
            @test MOI.is_valid(m, x[i])
            @test MOI.is_valid(m, x_int[i])
            @test MOI.is_valid(m, y[i])
            @test MOI.is_valid(m, y_bool[i])
            @test MOI.is_valid(m, z[i])
        end

        # Don't set names to check whether they are made unique before
        # generating the model.

        # Add constraints.
        c1 = MOI.add_constraint(
            m,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, 1], x), 0),
            MOI.EqualTo(2),
        )
        c2 = MOI.add_constraint(
            m,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, 1], x), 0),
            MOI.LessThan(2),
        )
        c3 = MOI.add_constraint(
            m,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, 1], y), 0),
            MOI.EqualTo(true),
        )
        c4 = MOI.add_constraint(
            m,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, 1], y), 0),
            MOI.EqualTo(0),
        )
        c5 = MOI.add_constraint(
            m,
            MOI.ScalarAffineFunction(
                MOI.ScalarAffineTerm.([1.0, 2.0], z),
                0.0,
            ),
            MOI.EqualTo(2.0),
        )
        c6 = MOI.add_constraint(
            m,
            MOI.ScalarAffineFunction(
                MOI.ScalarAffineTerm.([1.0, 2.0], z),
                0.0,
            ),
            MOI.LessThan(2.0),
        )
        c7 = MOI.add_constraint(
            m,
            MOI.ScalarAffineFunction(
                MOI.ScalarAffineTerm.([1.0, 2.0], z),
                0.0,
            ),
            CP.Strictly(MOI.LessThan(2.0)),
        )
        c8 = MOI.add_constraint(
            m,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, 1], y), 0),
            CP.DifferentFrom(true),
        )
        c9 = MOI.add_constraint(
            m,
            MOI.ScalarAffineFunction(
                MOI.ScalarAffineTerm.([1.0, 1.0], z),
                0.0,
            ),
            CP.DifferentFrom(1.0),
        )
        c10 = MOI.add_constraint(
            m,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, 1], x), 0),
            CP.DifferentFrom(1),
        )
        c11 = MOI.add_constraint(
            m,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, 1], y), 0),
            MOI.LessThan(1),
        )

        @test MOI.is_valid(m, c1)
        @test MOI.is_valid(m, c2)
        @test MOI.is_valid(m, c3)
        @test MOI.is_valid(m, c4)
        @test MOI.is_valid(m, c5)
        @test MOI.is_valid(m, c6)
        @test MOI.is_valid(m, c7)
        @test MOI.is_valid(m, c8)
        @test MOI.is_valid(m, c9)
        @test MOI.is_valid(m, c10)
        @test MOI.is_valid(m, c11)

        # Test some attributes for these constraints.
        @test length(MOI.get(m, MOI.ListOfConstraintTypesPresent())) == 11

        # Generate the FZN file.
        io = IOBuffer(truncate=true)
        write(io, m)
        fzn = String(take!(io))

        @test fzn == """var int: x1 :: output_var;
            var int: x2 :: output_var;
            var bool: x3 :: output_var;
            var bool: x4 :: output_var;
            var float: x5 :: output_var;
            var float: x6 :: output_var;

            constraint int_lin_eq([1, 1], [x1, x2], 2);
            constraint int_lin_le([1, 1], [x1, x2], 2);
            constraint bool_lin_eq([1, 1], [x3, x4], 1);
            constraint bool_lin_eq([1, 1], [x3, x4], 0);
            constraint float_lin_eq([1.0, 2.0], [x5, x6], 2.0);
            constraint float_lin_le([1.0, 2.0], [x5, x6], 2.0);
            constraint float_lin_lt([1.0, 2.0], [x5, x6], 2.0);
            constraint int_lin_ne([1, 1], [x3, x4], 1);
            constraint float_lin_ne([1.0, 1.0], [x5, x6], 1.0);
            constraint int_lin_ne([1, 1], [x1, x2], 1);
            constraint bool_lin_le([1, 1], [x3, x4], 1);

            solve satisfy;
            """

        # Test that the names have been correctly transformed.
        for v in [x..., y..., z...]
            vn = MOI.get(m, MOI.VariableName(), v)
            @test match(r"^x\d+$", vn) !== nothing
        end
    end

    @testset "Constraints: MOI.VariableIndex of integer in MOI.EqualTo / MOI.LessThan / CP.Strictly{MOI.LessThan} / CP.DifferentFrom" begin
        m = CP.FlatZinc.Model()
        @test MOI.is_empty(m)

        # Create variable.
        x, x_int = MOI.add_constrained_variable(m, MOI.Integer())
        y, y_bool = MOI.add_constrained_variable(m, MOI.ZeroOne())

        @test !MOI.is_empty(m)
        @test MOI.is_valid(m, x)
        @test MOI.is_valid(m, x_int)
        @test MOI.is_valid(m, y)
        @test MOI.is_valid(m, y_bool)

        # Don't set names to check whether they are made unique before
        # generating the model.

        # Add constraints.
        c1 = MOI.add_constraint(m, x, MOI.EqualTo(4))
        c2 = MOI.add_constraint(m, y, MOI.EqualTo(false))
        c3 = MOI.add_constraint(m, y, MOI.EqualTo(1)) # 1 will be cast to true.
        c4 = MOI.add_constraint(m, x, MOI.LessThan(2))
        c5 = MOI.add_constraint(m, y, MOI.LessThan(true))
        c6 = MOI.add_constraint(m, x, CP.Strictly(MOI.LessThan(2)))
        c7 = MOI.add_constraint(m, y, CP.Strictly(MOI.LessThan(true)))
        c8 = MOI.add_constraint(m, x, CP.DifferentFrom(2))
        c9 = MOI.add_constraint(m, y, CP.DifferentFrom(false))

        @test MOI.is_valid(m, c1)
        @test MOI.is_valid(m, c2)
        @test MOI.is_valid(m, c3)
        @test MOI.is_valid(m, c4)
        @test MOI.is_valid(m, c5)
        @test MOI.is_valid(m, c6)
        @test MOI.is_valid(m, c7)
        @test MOI.is_valid(m, c8)
        @test MOI.is_valid(m, c9)

        # Test some attributes for these constraints.
        @test length(MOI.get(m, MOI.ListOfConstraintTypesPresent())) == 10

        # Generate the FZN file.
        io = IOBuffer(truncate=true)
        write(io, m)
        fzn = String(take!(io))

        @test fzn == """var int: x1 :: output_var;
            var bool: x2 :: output_var;

            constraint int_eq(x1, 4);
            constraint bool_eq(x2, false);
            constraint bool_eq(x2, true);
            constraint int_le(x1, 2);
            constraint bool_le(x2, true);
            constraint int_lt(x1, 2);
            constraint bool_lt(x2, true);
            constraint int_ne(x1, 2);
            constraint int_ne(x2, false);

            solve satisfy;
            """

        # Test that the names have been correctly transformed.
        for v in [x, y]
            vn = MOI.get(m, MOI.VariableName(), v)
            @test match(r"^x\d+$", vn) !== nothing
        end
    end

    @testset "Constraints: MOI.VariableIndex of float in MOI.EqualTo / MOI.LessThan / CP.Strictly{MOI.LessThan} / CP.DifferentFrom" begin
        m = CP.FlatZinc.Model()
        @test MOI.is_empty(m)

        # Create variable.
        x = MOI.add_variable(m)

        @test !MOI.is_empty(m)
        @test MOI.is_valid(m, x)

        # Don't set names to check whether they are made unique before
        # generating the model.

        # Add constraints.
        c1 = MOI.add_constraint(m, x, MOI.EqualTo(5.0))
        c2 = MOI.add_constraint(m, x, MOI.LessThan(5.0))
        c3 = MOI.add_constraint(m, x, CP.Strictly(MOI.LessThan(5.0)))
        c4 = MOI.add_constraint(m, x, CP.DifferentFrom(5.0))

        @test MOI.is_valid(m, c1)
        @test MOI.is_valid(m, c2)
        @test MOI.is_valid(m, c3)
        @test MOI.is_valid(m, c4)

        # Test some attributes for these constraints.
        @test length(MOI.get(m, MOI.ListOfConstraintTypesPresent())) == 4

        # Generate the FZN file.
        io = IOBuffer(truncate=true)
        write(io, m)
        fzn = String(take!(io))

        @test fzn == """var float: x1 :: output_var;

            constraint float_eq(x1, 5.0);
            constraint float_le(x1, 5.0);
            constraint float_lt(x1, 5.0);
            constraint float_ne(x1, 5.0);

            solve satisfy;
            """

        # Test that the names have been correctly transformed.
        for v in [x]
            vn = MOI.get(m, MOI.VariableName(), v)
            @test match(r"^x\d+$", vn) !== nothing
        end
    end

    @testset "Constraints: CP.Reification{MOI.VectorOfVariables of integers in MOI.EqualTo / MOI.LessThan / CP.Strictly{MOI.LessThan} / CP.DifferentFrom}" begin
        m = CP.FlatZinc.Model()
        @test MOI.is_empty(m)

        # Create variable.
        x, x_bool = MOI.add_constrained_variable(m, MOI.ZeroOne())
        y, y_int = MOI.add_constrained_variable(m, MOI.Integer())

        @test !MOI.is_empty(m)
        @test MOI.is_valid(m, x)
        @test MOI.is_valid(m, x_bool)
        @test MOI.is_valid(m, y)
        @test MOI.is_valid(m, y_int)

        # Don't set names to check whether they are made unique before
        # generating the model.

        # Add constraints.
        c1 = MOI.add_constraint(m, [x, y], CP.Reification(MOI.EqualTo(2)))
        c2 = MOI.add_constraint(m, [x, y], CP.Reification(MOI.LessThan(2)))
        c3 = MOI.add_constraint(
            m,
            [x, y],
            CP.Reification(CP.Strictly(MOI.LessThan(2))),
        )
        c4 = MOI.add_constraint(m, [x, y], CP.Reification(CP.DifferentFrom(2)))

        @test MOI.is_valid(m, c1)
        @test MOI.is_valid(m, c2)
        @test MOI.is_valid(m, c3)
        @test MOI.is_valid(m, c4)

        # Test some attributes for these constraints.
        @test length(MOI.get(m, MOI.ListOfConstraintTypesPresent())) == 6

        # Generate the FZN file.
        io = IOBuffer(truncate=true)
        write(io, m)
        fzn = String(take!(io))

        @test fzn == """var bool: x1 :: output_var;
            var int: x2 :: output_var;

            constraint int_eq_reif(x2, 2, x1);
            constraint int_le_reif(x1, 2, x1);
            constraint int_lt_reif(x2, 2, x1);
            constraint int_ne_reif(x2, 2, x1);

            solve satisfy;
            """

        # Test that the names have been correctly transformed.
        for v in [x, y]
            vn = MOI.get(m, MOI.VariableName(), v)
            @test match(r"^x\d+$", vn) !== nothing
        end
    end

    @testset "Constraints: CP.Reification{MOI.VectorOfVariables of floats in MOI.EqualTo / MOI.LessThan / CP.Strictly{MOI.LessThan} / CP.DifferentFrom}" begin
        m = CP.FlatZinc.Model()
        @test MOI.is_empty(m)

        # Create variable.
        x, x_bool = MOI.add_constrained_variable(m, MOI.ZeroOne())
        y = MOI.add_variable(m)

        @test !MOI.is_empty(m)
        @test MOI.is_valid(m, x)
        @test MOI.is_valid(m, x_bool)
        @test MOI.is_valid(m, y)

        # Don't set names to check whether they are made unique before
        # generating the model.

        # Add constraints.
        c1 = MOI.add_constraint(m, [x, y], CP.Reification(MOI.EqualTo(2.0)))
        c2 = MOI.add_constraint(m, [x, y], CP.Reification(MOI.LessThan(2.0)))
        c3 = MOI.add_constraint(
            m,
            [x, y],
            CP.Reification(CP.Strictly(MOI.LessThan(2.0))),
        )
        c4 =
            MOI.add_constraint(m, [x, y], CP.Reification(CP.DifferentFrom(2.0)))

        @test MOI.is_valid(m, c1)
        @test MOI.is_valid(m, c2)
        @test MOI.is_valid(m, c3)
        @test MOI.is_valid(m, c4)

        # Test some attributes for these constraints.
        @test length(MOI.get(m, MOI.ListOfConstraintTypesPresent())) == 5

        # Generate the FZN file.
        io = IOBuffer(truncate=true)
        write(io, m)
        fzn = String(take!(io))

        @test fzn == """var bool: x1 :: output_var;
            var float: x2 :: output_var;

            constraint float_eq_reif(x2, 2.0, x1);
            constraint float_le_reif(x2, 2.0, x1);
            constraint float_lt_reif(x2, 2.0, x1);
            constraint float_ne_reif(x2, 2.0, x1);

            solve satisfy;
            """

        # Test that the names have been correctly transformed.
        for v in [x, y]
            vn = MOI.get(m, MOI.VariableName(), v)
            @test match(r"^x\d+$", vn) !== nothing
        end
    end

    @testset "Constraints: CP.Reification{MOI.VectorAffineFunction of integers in MOI.EqualTo / MOI.LessThan / CP.Strictly{MOI.LessThan} / CP.DifferentFrom}" begin
        m = CP.FlatZinc.Model()
        @test MOI.is_empty(m)

        # Create variable.
        w, w_bool = MOI.add_constrained_variable(m, MOI.ZeroOne())
        x, x_bool = MOI.add_constrained_variable(m, MOI.ZeroOne())
        y, y_int = MOI.add_constrained_variable(m, MOI.Integer())
        z, z_int = MOI.add_constrained_variable(m, MOI.Integer())

        @test !MOI.is_empty(m)
        @test MOI.is_valid(m, w)
        @test MOI.is_valid(m, w_bool)
        @test MOI.is_valid(m, x)
        @test MOI.is_valid(m, x_bool)
        @test MOI.is_valid(m, y)
        @test MOI.is_valid(m, y_int)
        @test MOI.is_valid(m, z)
        @test MOI.is_valid(m, z_int)

        # Don't set names to check whether they are made unique before
        # generating the model.

        # Add constraints.
        vaf = MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 2, 2, 2],
                MOI.ScalarAffineTerm.([1, 1, 1, 1], [w, x, y, z]),
            ),
            [0, 0],
        )

        c1 = MOI.add_constraint(m, vaf, CP.Reification(MOI.EqualTo(2)))
        c2 = MOI.add_constraint(m, vaf, CP.Reification(MOI.LessThan(2)))
        c3 = MOI.add_constraint(
            m,
            vaf,
            CP.Reification(CP.Strictly(MOI.LessThan(2))),
        )
        c4 = MOI.add_constraint(m, vaf, CP.Reification(CP.DifferentFrom(2)))

        @test MOI.is_valid(m, c1)
        @test MOI.is_valid(m, c2)
        @test MOI.is_valid(m, c3)
        @test MOI.is_valid(m, c4)

        # Test some attributes for these constraints.
        @test length(MOI.get(m, MOI.ListOfConstraintTypesPresent())) == 6

        # Generate the FZN file.
        io = IOBuffer(truncate=true)
        write(io, m)
        fzn = String(take!(io))

        @test fzn == """var bool: x1 :: output_var;
            var bool: x2 :: output_var;
            var int: x3 :: output_var;
            var int: x4 :: output_var;

            constraint int_lin_eq_reif([1, 1, 1], [x2, x3, x4], 2, x1);
            constraint int_lin_le_reif([1, 1, 1], [x2, x3, x4], 2, x1);
            constraint int_lin_lt_reif([1, 1, 1], [x2, x3, x4], 2, x1);
            constraint int_lin_ne_reif([1, 1, 1], [x2, x3, x4], 2, x1);

            solve satisfy;
            """

        # Test that the names have been correctly transformed.
        for v in [x, y]
            vn = MOI.get(m, MOI.VariableName(), v)
            @test match(r"^x\d+$", vn) !== nothing
        end
    end

    @testset "Constraints: CP.Reification{MOI.VectorAffineFunction of floats in MOI.EqualTo / MOI.LessThan / CP.Strictly{MOI.LessThan} / CP.DifferentFrom}" begin
        m = CP.FlatZinc.Model()
        @test MOI.is_empty(m)

        # Create variable.
        w, w_bool = MOI.add_constrained_variable(m, MOI.ZeroOne())
        x, x_bool = MOI.add_constrained_variable(m, MOI.ZeroOne())
        y, y_int = MOI.add_constrained_variable(m, MOI.Integer())
        z = MOI.add_variable(m)

        @test !MOI.is_empty(m)
        @test MOI.is_valid(m, w)
        @test MOI.is_valid(m, w_bool)
        @test MOI.is_valid(m, x)
        @test MOI.is_valid(m, x_bool)
        @test MOI.is_valid(m, y)
        @test MOI.is_valid(m, y_int)
        @test MOI.is_valid(m, z)

        # Don't set names to check whether they are made unique before
        # generating the model.

        # Add constraints.
        vaf = MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [1, 2, 2, 2],
                MOI.ScalarAffineTerm.([1.0, 1.0, 1.0, 1.0], [w, x, y, z]),
            ),
            [0.0, 0.0],
        )

        c1 = MOI.add_constraint(m, vaf, CP.Reification(MOI.EqualTo(2.0)))
        c2 = MOI.add_constraint(m, vaf, CP.Reification(MOI.LessThan(2.0)))
        c3 = MOI.add_constraint(
            m,
            vaf,
            CP.Reification(CP.Strictly(MOI.LessThan(2.0))),
        )
        c4 = MOI.add_constraint(m, vaf, CP.Reification(CP.DifferentFrom(2.0)))

        @test MOI.is_valid(m, c1)
        @test MOI.is_valid(m, c2)
        @test MOI.is_valid(m, c3)
        @test MOI.is_valid(m, c4)

        # Test some attributes for these constraints.
        @test length(MOI.get(m, MOI.ListOfConstraintTypesPresent())) == 6

        # Generate the FZN file.
        io = IOBuffer(truncate=true)
        write(io, m)
        fzn = String(take!(io))

        @test fzn == """var bool: x1 :: output_var;
            var bool: x2 :: output_var;
            var int: x3 :: output_var;
            var float: x4 :: output_var;

            constraint float_lin_eq_reif([1.0, 1.0, 1.0], [x2, x3, x4], 2.0, x1);
            constraint float_lin_le_reif([1.0, 1.0, 1.0], [x2, x3, x4], 2.0, x1);
            constraint float_lin_lt_reif([1.0, 1.0, 1.0], [x2, x3, x4], 2.0, x1);
            constraint float_lin_ne_reif([1.0, 1.0, 1.0], [x2, x3, x4], 2.0, x1);

            solve satisfy;
            """

        # Test that the names have been correctly transformed.
        for v in [x, y]
            vn = MOI.get(m, MOI.VariableName(), v)
            @test match(r"^x\d+$", vn) !== nothing
        end
    end

    @testset "Name rewriting" begin
        m = CP.FlatZinc.Model()
        @test MOI.is_empty(m)

        # Create variables.
        x, x_int = MOI.add_constrained_variable(m, MOI.Integer())
        y, y_bool =
            MOI.add_constrained_variables(m, [MOI.ZeroOne() for _ in 1:5])

        # Set names. The name for x in invalid for FlatZinc, but not
        # the ones for y.
        MOI.set(m, MOI.VariableName(), x, "_x")
        for i in 1:5
            MOI.set(m, MOI.VariableName(), y[i], "y_$(i)")
        end

        @test MOI.get(m, MOI.VariableName(), x) == "_x"
        for i in 1:5
            @test MOI.get(m, MOI.VariableName(), y[i]) == "y_$(i)"
        end

        # Generate the FZN file.
        io = IOBuffer(truncate=true)
        write(io, m)
        fzn = String(take!(io))

        @test fzn == """var int: x_x :: output_var;
            var bool: y_1 :: output_var;
            var bool: y_2 :: output_var;
            var bool: y_3 :: output_var;
            var bool: y_4 :: output_var;
            var bool: y_5 :: output_var;

            solve satisfy;
            """

        # Test that the names have been correctly transformed.
        @test MOI.get(m, MOI.VariableName(), x) == "x_x"
        for i in 1:5
            @test MOI.get(m, MOI.VariableName(), y[i]) == "y_$(i)"
        end
    end

    @testset "Maximising" begin
        m = CP.FlatZinc.Model()
        @test MOI.is_empty(m)

        x, x_int = MOI.add_constrained_variable(m, MOI.Integer())

        @test !MOI.is_empty(m)
        @test MOI.is_valid(m, x)
        @test MOI.is_valid(m, x_int)

        MOI.set(
            m,
            MOI.ObjectiveFunction{MOI.VariableIndex}(),
            x,
        )

        MOI.set(m, MOI.ObjectiveSense(), MOI.MAX_SENSE)

        io = IOBuffer(truncate=true)
        write(io, m)
        fzn = String(take!(io))

        @test fzn == """var int: x1 :: output_var;

        solve maximize x1;
        """
    end

    @testset "Minimising" begin
        m = CP.FlatZinc.Model()
        @test MOI.is_empty(m)

        x, x_int = MOI.add_constrained_variable(m, MOI.Integer())

        @test !MOI.is_empty(m)
        @test MOI.is_valid(m, x)
        @test MOI.is_valid(m, x_int)

        MOI.set(
            m,
            MOI.ObjectiveFunction{MOI.VariableIndex}(),
            x,
        )

        MOI.set(m, MOI.ObjectiveSense(), MOI.MIN_SENSE)

        io = IOBuffer(truncate=true)
        write(io, m)
        fzn = String(take!(io))

        @test fzn == """var int: x1 :: output_var;

        solve minimize x1;
        """
    end
end


text:
using LibSerialPort

include("SerialDevice.jl")


text:
using Tethys
using Printf

MAX_ORDER=2
diag=Tethys.Diag(0,1.0,MAX_ORDER,[0.1,0.2,1.0 , 0.3,0.4,1.0], [0,0,0,0,0,0])

GF=Tethys.BareExpansion(diag)

println("Splish splash: ",diag," ==♓> ", GF)

GF=Tethys.Monte!(diag)
println("Monte: ", GF)

fout=open("Monte.dat","w")
for i in 1:100_000
    GF=Tethys.Monte!(diag)
    @printf(fout,"%d %f \n",i,GF)
end
close(fout)

println("Final diag: $(diag) GF: $(GF)")


text:
# Tests for General Direct Approach Cost DEA Models
@testset "CostGDADEAModel" begin

    # Test using Book data
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6 8];
    Y = [1; 1; 1; 1; 1; 1; 1; 1];
    W = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1];

    # Cost GDA
    costgda = deacostgda(X, Y, W, :ERG)

    @test typeof(costgda) == CostGDADEAModel

    @test nobs(costgda) == 8
    @test ninputs(costgda) == 2
    @test noutputs(costgda) == 1
    @test ismonetary(costgda) == false

    @test efficiency(costgda) == efficiency(costgda, :Economic)
    @test efficiency(costgda, :Economic)   ≈ [0.0; 0.5; 0.5; 5/12; 0.6; 0.25; 0.525; 0.532609] atol = 1e-5
    @test efficiency(costgda, :Technical)  ≈ [0.0; 0.0; 0.0; 5/12; 0.6; 1/6; 0.35; 7/16]
    @test efficiency(costgda, :Allocative) ≈ [0.0; 0.5; 0.5; 0.0; 0.0; 1/12; 0.175; 0.095109] atol = 1e-5
    @test normfactor(costgda) ≈ [4.0; 2.0; 2.0; 7.2; 10.0; 12.0; 40/7; 368/35]

    # With monetary = true
    costgdamon = deacostgda(X, Y, W, :ERG, monetary = true)

    @test ismonetary(costgdamon) == true

    @test efficiency(costgdamon) ≈ efficiency(costgda) .* normfactor(costgda)
    @test efficiency(costgdamon, :Technical) ≈ efficiency(costgda, :Technical) .* normfactor(costgda)
    @test efficiency(costgdamon, :Allocative) ≈ efficiency(costgda, :Allocative) .* normfactor(costgda)

    # With CRS (Same result in this example)
    costgdacrs = deacostgda(X, Y, W, :ERG, rts = :CRS)

    @test efficiency(costgdacrs, :Economic)   ≈ [0.0; 0.5; 0.5; 5/12; 0.6; 0.25; 0.525; 0.532609] atol = 1e-5
    @test efficiency(costgdacrs, :Technical)  ≈ [0.0; 0.0; 0.0; 5/12; 0.6; 1/6; 0.35; 7/16]
    @test efficiency(costgdacrs, :Allocative) ≈ [0.0; 0.5; 0.5; 0.0; 0.0; 1/12; 0.175; 0.095109] atol = 1e-5

    # Print
    show(IOBuffer(), costgda)

    # Test errors
    @test_throws DimensionMismatch deacostgda([1; 2 ; 3], [4 ; 5], [1; 1; 1], :ERG) #  Different number of observations
    @test_throws DimensionMismatch deacostgda([1; 2; 3], [4; 5; 6], [1; 2; 3; 4], :ERG) # Different number of observation in prices
    @test_throws DimensionMismatch deacostgda([1 1; 2 2; 3 3 ], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], :ERG) # Different number of input prices and inputs
    @test_throws ArgumentError deacostgda([1; 2; 3], [4; 5; 6], [1; 2; 3], :Error) # Invalid associated efficiency measure
    @test_throws ArgumentError deacostgda([1; 2; 3], [4; 5; 6], [1; 2; 3], :ERG, rts = :Error) # Invalid returns to scale

end


text:
"""This module defines classes of the package."""
# Author: Henri Gérard <hgerard.proy@gmail.com>
# License: MIT

abstract type RegressionModel end

# Define an object to
type LinearRegression <: RegressionModel
    function LinearRegression()
        return new()
    end
end

type LogisticRegression <: RegressionModel
    function LogisticRegression()
        return new()
    end
end

abstract type ParallelParam end

type Sequential <: ParallelParam
    function Sequential()
        return new()
    end
end

type Parallel <: ParallelParam
    function Parallel()
        return new()
    end
end

abstract type AlgoParams end

type OptParams <: AlgoParams
    itmax::Int
    stability::Float64
    learning_rate::Float64
    verbosity::Int
    function OptParams(itmax, stability, learning_rate; verbosity = 0)
        return new(itmax, stability, learning_rate, verbosity)
    end
end

type ProjParams <: AlgoParams
    ITER_MAX::Int
    precision::Float64
    sample::Int
    para_proj::ParallelParam
    para_inter::ParallelParam
    function ProjParams(ITER_MAX, precision, sample; para_proj=Sequential(), para_inter=Sequential())
        return new(ITER_MAX, precision, sample, para_proj, para_inter)
    end
end

abstract type AmbiguitySet end

type DivergenceSet <: AmbiguitySet
    function DivergenceSet()
        return new()
    end
end

type WassersteinSet <: AmbiguitySet
    function WassersteinSet()
        return new()
    end
end

abstract type GeneralConstraint end

type PositiveConstraint <: GeneralConstraint
    function PositiveConstraint()
        return new()
    end
end

type DROConstraint <: GeneralConstraint
    function DROConstraint()
        return new()
    end
end

type EntropicConstraint <: GeneralConstraint
    function EntropicConstraint()
        return new()
    end
end

abstract type DivergenceConstraint<:GeneralConstraint end

type KLConstraint <: DivergenceConstraint
    function KLConstraint()
        return new()
    end
end

type RobustModel
    descent_direction::Array{Float64,1}
    I0::UnitRange{Int64}
    name::GeneralConstraint
    regressionModel::RegressionModel
    function RobustModel(N, nb_features, ϵ, ambiguity, regressionModel)
        if ambiguity == "KLdivergence"
            descent_direction = [ϵ; 1; zeros(1:nb_features); ones(N)/N]
            I0 = 1:N
            dim = N+1
            name = KLConstraint()
        elseif ambiguity == "wasserstein"
            descent_direction = [ϵ; zeros(1:nb_features); ones(N)/N]
            I0 = 1:N^2
            dim = N^2+1
            name = DROConstraint()
        elseif ambiguity == "entropic"
            descent_direction = [1; zeros(1:nb_features); ones(N)/N]
            I0 = 1:N
            dim = N+1
            name = EntropicConstraint()
        end
        return new(descent_direction, I0, name, regressionModel)
    end
end

max_stars_repo_name: Henri-Gerard/robox

max_issues_repo_name: Henri-Gerard/robox

max_forks_repo_name: Henri-Gerard/robox


text:
using Test, PolynomialBases
import SymPy

x, α, β, N = SymPy.symbols("x, alpha, beta, N")

@test 0 == SymPy.simplify( hahn(x, 0, α, β, N) - 1 )
# https://www.wolframalpha.com/input/?i=HypergeometricPFQ%5B+%7B-1,+1%2Balpha%2Bbeta%2B1,+-x%7D,+%7Balpha%2B1,+-N%7D,+1%5D
@test 0 == SymPy.simplify( hahn(x, 1, α, β, N) - ( 1 - (α+β+2)*x/(N*(α+1)) ) )
# https://www.wolframalpha.com/input/?i=HypergeometricPFQ%5B+%7B-2,+2%2Balpha%2Bbeta%2B1,+-x%7D,+%7Balpha%2B1,+-N%7D,+1%5D
@test 0 == SymPy.simplify( hahn(x, 2, α, β, N) - (
        (1-x)*x*(α+β+4)*(α+β+3) / ((α+1)*(α+2)*(1-N)*N)
        - 2x*(α+β+3)/(N*(α+1)) + 1
    ) )
# https://www.wolframalpha.com/input/?i=HypergeometricPFQ%5B+%7B-3,+3%2Balpha%2Bbeta%2B1,+-x%7D,+%7Balpha%2B1,+-N%7D,+1%5D
@test 0 == SymPy.simplify( hahn(x, 3, α, β, N) - (
        3*(1-x)*x*(α+β+5)*(α+β+4) / ( (α+1)*(α+2)*N*(1-N) )
        - (1-x)*(2-x)*x*(α+β+5)*(α+β+6)*(α+β+4) / ( (α+1)*(α+2)*(α+3)*N*(1-N)*(2-N) )
        - 3x*(α+β+4) / ( N*(α+1) ) + 1
    ) )
# https://www.wolframalpha.com/input/?i=HypergeometricPFQ%5B+%7B-4,+4%2Balpha%2Bbeta%2B1,+-x%7D,+%7Balpha%2B1,+-N%7D,+1%5D
@test 0 == SymPy.simplify( hahn(x, 4, α, β, N) - (
        6*(1-x)*x*(α+β+6)*(α+β+5) / ( (α+1)*(α+2)*N*(1-N) )
        - 4*(1-x)*(2-x)*x*(α+β+6)*(α+β+7)*(α+β+5) / ( (α+1)*(α+2)*(α+3)*N*(1-N)*(2-N) )
        + (1-x)*(2-x)*(3-x)*x*(α+β+6)*(α+β+7)*(α+β+8)*(α+β+5) / ( (α+1)*(α+2)*(α+3)*(α+4)*N*(1-N)*(2-N)*(3-N) )
        - 4x*(α+β+5) / ( N*(α+1) )
        + 1
    ) )

@test_throws ArgumentError hahn(x, 6, α, β, 5)
@inferred hahn(x, 4, α, β, 5)
@inferred hahn(10., 4, -0.5, -0.9, 5)


text:
import Compat.LinearAlgebra: Givens, AbstractRotation

### These three mul! should be in Base, but for the time being they do not add methods to Base.mul!; they add methods to the internal mul!.

function mul!(G::Givens{T}, A::AbstractVecOrMat) where T<:Real
    m, n = size(A, 1), size(A, 2)
    if G.i2 > m
        throw(DimensionMismatch("column indices for rotation are outside the matrix"))
    end
    @inbounds @simd for i = 1:n
        a1, a2 = A[G.i1,i], A[G.i2,i]
        A[G.i1,i] =  G.c*a1 + G.s*a2
        A[G.i2,i] = -G.s*a1 + G.c*a2
    end
    return A
end

function mul!(A::AbstractMatrix, G::Givens)
    m, n = size(A, 1), size(A, 2)
    if G.i2 > n
        throw(DimensionMismatch("column indices for rotation are outside the matrix"))
    end
    @inbounds @simd for i = 1:m
        a1, a2 = A[i,G.i1], A[i,G.i2]
        A[i,G.i1] = G.c*a1 - conj(G.s)*a2
        A[i,G.i2] = G.s*a1 +       G.c*a2
    end
    return A
end

function mul!(A::AbstractMatrix, G::Givens{T}) where T<:Real
    m, n = size(A, 1), size(A, 2)
    if G.i2 > n
        throw(DimensionMismatch("column indices for rotation are outside the matrix"))
    end
    @inbounds @simd for i = 1:m
        a1, a2 = A[i,G.i1], A[i,G.i2]
        A[i,G.i1] = G.c*a1 - G.s*a2
        A[i,G.i2] = G.s*a1 + G.c*a2
    end
    return A
end

struct Pnmp2toPlm{T} <: AbstractRotation{T}
    rotations::Vector{Givens{T}}
end

function Pnmp2toPlm(::Type{T}, n::Int, m::Int) where T
    G = Vector{Givens{T}}(undef, n)
    @inbounds for ℓ = 1:n
        c = sqrt(T((2m+2)*(2ℓ+2m+3))/T((ℓ+2m+2)*(ℓ+2m+3)))
        s = sqrt(T(ℓ*(ℓ+1))/T((ℓ+2m+2)*(ℓ+2m+3)))
        G[n+1-ℓ] = Givens(ℓ, ℓ+2, c, s)
    end
    Pnmp2toPlm(G)
end

@inline length(P::Pnmp2toPlm) = length(P.rotations)
@inline getindex(P::Pnmp2toPlm, i::Int) = P.rotations[i]

function LAmul!(C::Pnmp2toPlm, A::AbstractVecOrMat)
    @inbounds for i = 1:length(C)
        mul!(C.rotations[i], A)
    end
    A
end

function LAmul!(A::AbstractMatrix, C::Pnmp2toPlm)
    @inbounds for i = length(C):-1:1
        mul!(A, C.rotations[i])
    end
    A
end

struct RotationPlan{T} <: AbstractRotation{T}
    layers::Vector{Pnmp2toPlm{T}}
    snm::Vector{T}
    cnm::Vector{T}
end

function RotationPlan(::Type{T}, n::Int) where T
    layers = Vector{Pnmp2toPlm{T}}(undef, n-1)
    @inbounds for m = 0:n-2
        layers[m+1] = Pnmp2toPlm(T, n-1-m, m)
    end
    n = n+1
    snm = zeros(T, (n*(n+1))÷2)
    cnm = zeros(T, (n*(n+1))÷2)
    @inbounds for l = 0:n-1
        for m = 0:n-l-1
            nums = T((l+1)*(l+2))
            numc = T((2*m+2)*(2*l+2*m+5))
            den = T((l+2*m+3)*(l+2*m+4))
            snm[l+(m*(2*n+1-m))÷2+1] = sqrt(nums/den)
            cnm[l+(m*(2*n+1-m))÷2+1] = sqrt(numc/den)
        end
    end
    RotationPlan(layers, snm, cnm)
end

if VERSION ≥ v"0.7-"
    adjoint(P::RotationPlan) = Adjoint(P)
    transpose(P::RotationPlan) = Transpose(P)
end

function LAmul!(P::RotationPlan, A::AbstractMatrix)
    N, M = size(A)
    snm = P.snm
    cnm = P.cnm
    if isodd(M)
        m = M÷2
        @inbounds for j = m:-2:2
            for l = N-j:-1:1
                s = snm[l+(j-2)*(2*N+3-j)÷2]
                c = cnm[l+(j-2)*(2*N+3-j)÷2]
                a1 = A[l+N*(2*m-1)]
                a2 = A[l+2+N*(2*m-1)]
                a3 = A[l+N*(2*m)]
                a4 = A[l+2+N*(2*m)]
                A[l+N*(2*m-1)] = c*a1 + s*a2
                A[l+2+N*(2*m-1)] = c*a2 - s*a1
                A[l+N*(2*m)] = c*a3 + s*a4
                A[l+2+N*(2*m)] = c*a4 - s*a3
            end
        end
    else
        m = M÷2
        @inbounds for j = m:-2:2
            for l = N-j:-1:1
                s = snm[l+(j-2)*(2*N+3-j)÷2]
                c = cnm[l+(j-2)*(2*N+3-j)÷2]
                a1 = A[l+N*(2*m-1)]
                a2 = A[l+2+N*(2*m-1)]
                A[l+N*(2*m-1)] = c*a1 + s*a2
                A[l+2+N*(2*m-1)] = c*a2 - s*a1
            end
        end
    end
    @stepthreads for m = M÷2-1:-1:2
        @inbounds for j = m:-2:2
            for l = N-j:-1:1
                s = snm[l+(j-2)*(2*N+3-j)÷2]
                c = cnm[l+(j-2)*(2*N+3-j)÷2]
                a1 = A[l+N*(2*m-1)]
                a2 = A[l+2+N*(2*m-1)]
                a3 = A[l+N*(2*m)]
                a4 = A[l+2+N*(2*m)]
                A[l+N*(2*m-1)] = c*a1 + s*a2
                A[l+2+N*(2*m-1)] = c*a2 - s*a1
                A[l+N*(2*m)] = c*a3 + s*a4
                A[l+2+N*(2*m)] = c*a4 - s*a3
            end
        end
    end
    A
end

if VERSION < v"0.7-"
    function Base.At_mul_B!(P::RotationPlan, A::AbstractMatrix)
        N, M = size(A)
        snm = P.snm
        cnm = P.cnm
        if isodd(M)
            m = M÷2
            @inbounds for j = reverse(m:-2:2)
                for l = 1:N-j
                    s = snm[l+(j-2)*(2*N+3-j)÷2]
                    c = cnm[l+(j-2)*(2*N+3-j)÷2]
                    a1 = A[l+N*(2*m-1)]
                    a2 = A[l+2+N*(2*m-1)]
                    a3 = A[l+N*(2*m)]
                    a4 = A[l+2+N*(2*m)]
                    A[l+N*(2*m-1)] = c*a1 - s*a2
                    A[l+2+N*(2*m-1)] = c*a2 + s*a1
                    A[l+N*(2*m)] = c*a3 - s*a4
                    A[l+2+N*(2*m)] = c*a4 + s*a3
                end
            end
        else
            m = M÷2
            @inbounds for j = reverse(m:-2:2)
                for l = 1:N-j
                    s = snm[l+(j-2)*(2*N+3-j)÷2]
                    c = cnm[l+(j-2)*(2*N+3-j)÷2]
                    a1 = A[l+N*(2*m-1)]
                    a2 = A[l+2+N*(2*m-1)]
                    A[l+N*(2*m-1)] = c*a1 - s*a2
                    A[l+2+N*(2*m-1)] = c*a2 + s*a1
                end
            end
        end
        @stepthreads for m = M÷2-1:-1:2
            @inbounds for j = reverse(m:-2:2)
                for l = 1:N-j
                    s = snm[l+(j-2)*(2*N+3-j)÷2]
                    c = cnm[l+(j-2)*(2*N+3-j)÷2]
                    a1 = A[l+N*(2*m-1)]
                    a2 = A[l+2+N*(2*m-1)]
                    a3 = A[l+N*(2*m)]
                    a4 = A[l+2+N*(2*m)]
                    A[l+N*(2*m-1)] = c*a1 - s*a2
                    A[l+2+N*(2*m-1)] = c*a2 + s*a1
                    A[l+N*(2*m)] = c*a3 - s*a4
                    A[l+2+N*(2*m)] = c*a4 + s*a3
                end
            end
        end
        A
    end
    Base.Ac_mul_B!(P::RotationPlan, A::AbstractMatrix) = At_mul_B!(P, A)
else
    function LinearAlgebra.lmul!(Pt::Transpose{T,<:RotationPlan}, A::AbstractMatrix) where T
        P = parent(Pt)
        N, M = size(A)
        snm = P.snm
        cnm = P.cnm
        if isodd(M)
            m = M÷2
            @inbounds for j = reverse(m:-2:2)
                for l = 1:N-j
                    s = snm[l+(j-2)*(2*N+3-j)÷2]
                    c = cnm[l+(j-2)*(2*N+3-j)÷2]
                    a1 = A[l+N*(2*m-1)]
                    a2 = A[l+2+N*(2*m-1)]
                    a3 = A[l+N*(2*m)]
                    a4 = A[l+2+N*(2*m)]
                    A[l+N*(2*m-1)] = c*a1 - s*a2
                    A[l+2+N*(2*m-1)] = c*a2 + s*a1
                    A[l+N*(2*m)] = c*a3 - s*a4
                    A[l+2+N*(2*m)] = c*a4 + s*a3
                end
            end
        else
            m = M÷2
            @inbounds for j = reverse(m:-2:2)
                for l = 1:N-j
                    s = snm[l+(j-2)*(2*N+3-j)÷2]
                    c = cnm[l+(j-2)*(2*N+3-j)÷2]
                    a1 = A[l+N*(2*m-1)]
                    a2 = A[l+2+N*(2*m-1)]
                    A[l+N*(2*m-1)] = c*a1 - s*a2
                    A[l+2+N*(2*m-1)] = c*a2 + s*a1
                end
            end
        end
        @stepthreads for m = M÷2-1:-1:2
            @inbounds for j = reverse(m:-2:2)
                for l = 1:N-j
                    s = snm[l+(j-2)*(2*N+3-j)÷2]
                    c = cnm[l+(j-2)*(2*N+3-j)÷2]
                    a1 = A[l+N*(2*m-1)]
                    a2 = A[l+2+N*(2*m-1)]
                    a3 = A[l+N*(2*m)]
                    a4 = A[l+2+N*(2*m)]
                    A[l+N*(2*m-1)] = c*a1 - s*a2
                    A[l+2+N*(2*m-1)] = c*a2 + s*a1
                    A[l+N*(2*m)] = c*a3 - s*a4
                    A[l+2+N*(2*m)] = c*a4 + s*a3
                end
            end
        end
        A
    end
    LinearAlgebra.lmul!(Pc::Adjoint{T,<:RotationPlan}, A::AbstractMatrix) where T = lmul!(parent(Pc)', A)

end

#=
function LAmul!(P::RotationPlan, A::AbstractMatrix)
    M, N = size(A)
    @inbounds for m = N÷2-2:-1:0
        layer = P.layers[m+1]
        for ℓ = 2*(m+2):4:N
            @simd for i = 1:length(layer)
                G = layer[i]
                a1, a2 = A[G.i1,ℓ], A[G.i2,ℓ]
                A[G.i1,ℓ] = G.c*a1 + G.s*a2
                A[G.i2,ℓ] = G.c*a2 - G.s*a1
                a1, a2 = A[G.i1,ℓ+1], A[G.i2,ℓ+1]
                A[G.i1,ℓ+1] = G.c*a1 + G.s*a2
                A[G.i2,ℓ+1] = G.c*a2 - G.s*a1
            end
        end
    end
    A
end

function Base.At_mul_B!(P::RotationPlan, A::AbstractMatrix)
    M, N = size(A)
    @inbounds for m = 0:N÷2-2
        layer = P.layers[m+1]
        for ℓ = 2*(m+2):4:N
            @simd for i = length(layer):-1:1
                G = layer[i]
                a1, a2 = A[G.i1,ℓ], A[G.i2,ℓ]
                A[G.i1,ℓ] = G.c*a1 - G.s*a2
                A[G.i2,ℓ] = G.s*a1 + G.c*a2
                a1, a2 = A[G.i1,ℓ+1], A[G.i2,ℓ+1]
                A[G.i1,ℓ+1] = G.c*a1 - G.s*a2
                A[G.i2,ℓ+1] = G.s*a1 + G.c*a2
            end
        end
    end
    A
end
=#

struct SlowSphericalHarmonicPlan{T} <: SphericalHarmonicPlan{T}
    RP::RotationPlan{T}
    p1::NormalizedLegendreToChebyshevPlan{T}
    p2::NormalizedLegendre1ToChebyshev2Plan{T}
    p1inv::ChebyshevToNormalizedLegendrePlan{T}
    p2inv::Chebyshev2ToNormalizedLegendre1Plan{T}
    B::Matrix{T}
end

function SlowSphericalHarmonicPlan(A::Matrix{T}) where T
    M, N = size(A)
    RP = RotationPlan(T, M-1)
    p1 = plan_normleg2cheb(A)
    p2 = plan_normleg12cheb2(A)
    p1inv = plan_cheb2normleg(A)
    p2inv = plan_cheb22normleg1(A)
    B = zero(A)
    SlowSphericalHarmonicPlan(RP, p1, p2, p1inv, p2inv, B)
end

function LAmul!(Y::Matrix, SP::SlowSphericalHarmonicPlan, X::Matrix)
    RP, p1, p2, B = SP.RP, SP.p1, SP.p2, SP.B
    copyto!(B, X)
    mul!(RP, B)
    M, N = size(X)
    mul_col_J!!(Y, p1, B, 1)
    @stepthreads for J = 2:4:N
        mul_col_J!!(Y, p2, B, J)
        J < N && mul_col_J!!(Y, p2, B, J+1)
    end
    @stepthreads for J = 4:4:N
        mul_col_J!!(Y, p1, B, J)
        J < N && mul_col_J!!(Y, p1, B, J+1)
    end
    Y
end

if  VERSION < v"0.7-"
    function Base.At_mul_B!(Y::Matrix, SP::SlowSphericalHarmonicPlan, X::Matrix)
        RP, p1inv, p2inv, B = SP.RP, SP.p1inv, SP.p2inv, SP.B
        copyto!(B, X)
        M, N = size(X)
        mul_col_J!!(Y, p1inv, B, 1)
        @stepthreads for J = 2:4:N
            mul_col_J!!(Y, p2inv, B, J)
            J < N && mul_col_J!!(Y, p2inv, B, J+1)
        end
        @stepthreads for J = 4:4:N
            mul_col_J!!(Y, p1inv, B, J)
            J < N && mul_col_J!!(Y, p1inv, B, J+1)
        end
        sph_zero_spurious_modes!(At_mul_B!(RP, Y))
    end

    Base.Ac_mul_B!(Y::Matrix, SP::SlowSphericalHarmonicPlan, X::Matrix) = At_mul_B!(Y, SP, X)
else
    transpose(P::SlowSphericalHarmonicPlan) = Transpose(P)
    adjoint(P::SlowSphericalHarmonicPlan) = Adjoint(P)

    function LinearAlgebra.mul!(Y::Matrix, SPt::Transpose{T,<:SlowSphericalHarmonicPlan}, X::Matrix) where T
        SP = parent(SPt)
        RP, p1inv, p2inv, B = SP.RP, SP.p1inv, SP.p2inv, SP.B
        copyto!(B, X)
        M, N = size(X)
        mul_col_J!!(Y, p1inv, B, 1)
        @stepthreads for J = 2:4:N
            mul_col_J!!(Y, p2inv, B, J)
            J < N && mul_col_J!!(Y, p2inv, B, J+1)
        end
        @stepthreads for J = 4:4:N
            mul_col_J!!(Y, p1inv, B, J)
            J < N && mul_col_J!!(Y, p1inv, B, J+1)
        end
        sph_zero_spurious_modes!(lmul!(transpose(RP), Y))
    end

    LinearAlgebra.mul!(Y::Matrix, SPc::Adjoint{T,<:SlowSphericalHarmonicPlan}, X::Matrix) where T =
        mul!(Y, transpose(parent(SPc)), X)
end


text:
function AbstractMultipleLinearRegression()
    return AbstractMultipleLinearRegression(())
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function estimate_error_variance(obj::AbstractMultipleLinearRegression)
    return jcall(obj, "estimateErrorVariance", jdouble, ())
end

function estimate_regressand_variance(obj::AbstractMultipleLinearRegression)
    return jcall(obj, "estimateRegressandVariance", jdouble, ())
end

function estimate_regression_parameters(obj::AbstractMultipleLinearRegression)
    return jcall(obj, "estimateRegressionParameters", Vector{jdouble}, ())
end

function estimate_regression_parameters_standard_errors(obj::AbstractMultipleLinearRegression)
    return jcall(obj, "estimateRegressionParametersStandardErrors", Vector{jdouble}, ())
end

function estimate_regression_parameters_variance(obj::AbstractMultipleLinearRegression)
    return jcall(obj, "estimateRegressionParametersVariance", Vector{Vector{jdouble}}, ())
end

function estimate_regression_standard_error(obj::AbstractMultipleLinearRegression)
    return jcall(obj, "estimateRegressionStandardError", jdouble, ())
end

function estimate_residuals(obj::AbstractMultipleLinearRegression)
    return jcall(obj, "estimateResiduals", Vector{jdouble}, ())
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function is_no_intercept(obj::AbstractMultipleLinearRegression)
    return jcall(obj, "isNoIntercept", jboolean, ())
end

function new_sample_data(obj::AbstractMultipleLinearRegression, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(obj, "newSampleData", void, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function set_no_intercept(obj::AbstractMultipleLinearRegression, arg0::jboolean)
    return jcall(obj, "setNoIntercept", void, (jboolean,), arg0)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
module LB

export Playground, FluidKernel, SingleRelaxationTime
export NashZeroOrderPressure, FluidStreaming, HalfWayBounceBack, ConstantVelocityIOlet
export ConstantPopulationIOlet, density, momentum, velocity, ParabolicVelocityIOlet

using ..AbstractLattice
using ..Simulation
using ..Playground: Feature, NOTHING
using ..Indices: GridCoords, index, Indexing, gridcoords
import ..Playground.initialize
import ..Indices.neighbor_index
import Base: length, ndims

# Base type for all kernally stuff
abstract Kernel
# Local kernel for each lattice type
abstract LocalKernel <: Kernel
# Base type for all collision kernels
abstract Collision <: Kernel
# Base type for all collision kernels
abstract Streaming <: Kernel
# Base type for all initilizers
abstract Initializer <: Kernel
# Holds local quantities
# - density
# - momentum
# - velocity
# - feq
abstract AbstractLocalQuantities

type FluidKernel <: LocalKernel
    # Collision kernel type and data
    collision::Collision
    # Describes how streaming takes place
    streamers :: Dict{Feature, Streaming}
end
type NullKernel <: LocalKernel; end

include("lattice.jl")
include("thermodynamics.jl")
include("LocalQuantities.jl")
include("collision.jl")
include("streaming.jl")
include("iolet.jl")

local_kernel(kernel::NullKernel, args...; kargs...) = nothing
function local_kernel(kernel::LocalKernel, sim::Simulation, site::Integer)
    quantities = LocalQuantities(sim.lattice)
    fᵢ = similar(sim.lattice.weights)
    local_kernel(kernel, sim, site, quantities, fᵢ)
end

function local_kernel(
    kernel::LocalKernel, sim::Simulation, site::Integer, quantities::LocalQuantities,
    fᵢ::DenseVector)
    @assert site > 0 && site < length(sim.indexing)
    @inbounds fᵢ[:] = sim.populations[:, site]
    @inbounds LocalQuantities!(quantities, fᵢ, sim.lattice)
    @inbounds collision!(sim.populations, kernel.collision, quantities.feq, site)
    for direction in 1:length(sim.lattice.weights)
        @inbounds const to = neighbor_index(sim, site, direction)
        const link = to == 0 ? NOTHING: sim.playground[to]
        const streamer = get(kernel.streamers, link, NULLSTREAMER)
        streaming(streamer, quantities, sim, site, to, direction)
    end
end

type Homogeneous{T <: AbstractFloat} <: Initializer
    density::T
    momentum::DenseVector{T}
end
type NullInitializer <: Initializer; end
const NULLINITIALIZER = NullInitializer()

function initialize(init::Homogeneous, sim::Simulation, index::Integer)
    sim.populations[:, index] = equilibrium(sim.lattice, init.density, init.momentum)
end
initialize(init::Homogeneous, sim::Simulation, coords::GridCoords) =
    initialize(init, sim, index(sim.indexing, coords))
initialize(::NullInitializer, sim::Simulation, index::Integer) = nothing

# Loop over all fluid sites and initialize populations
function initialize(sim::Simulation)
    for (site, feature) in enumerate(sim.playground)
        if haskey(sim.initializers, feature)
            initialize(sim.initializers[feature], sim, site)
        end
    end
end
# Loop over all fluid sites and perform calculations
function local_kernel(sim::Simulation)
    quants = LocalQuantities(sim.lattice)
    fᵢ = similar(sim.lattice.weights)
    for (site, feature) in enumerate(sim.playground)
        if haskey(sim.kernels, feature)
            local_kernel(sim.kernels[feature], sim, site, quants, fᵢ)
        end
    end
end
end


text:
function _trim(x::Array{T})::Tuple{Array{Int64},Array{Int64}} where {T<:AbstractFloat}
    ids = sortperm(x)
    distances = [x[ids[i+1]] - x[ids[i]] for i=1:length(x)-1]

    if length(x) < 3
        return Array{Int64, 1}(0), collect(1:length(x))
    end

    density = zeros(length(x))
    for i in eachindex(density)
        if i == 1
            density[i] = 0.25 / distances[1]
        elseif i == length(density)
            density[i] = 0.25 / distances[end]
        else
            density[i] = 1.0 / (distances[i-1] + distances[i])
        end
    end
    #normalize density
    density /= sum(density)

    #remove tails -> trim to 1σ
    cutoff = (1 - 0.6827) * 0.5

    prob = 0.0

    start = 1
    stop = length(x)

    for i in eachindex(x)
        prob += density[i]
        if prob >= cutoff
            start = max(1, i - 1)
            break
        end
    end

    j = length(x)
    prob = 0.0
    for i in eachindex(x)
        prob += density[j]
        if prob >= cutoff
            stop = min(length(x), j + 1)
            break
        end
        j -= 1
    end

    sort!(ids[[1:start..., stop:length(x)...]]), sort!(ids[[start+1:stop-1...]])
end

function trim(x::Array)
    deleteat!(x, _trim(x)[1])
    x
end

function trim(
    res::IntermediateResults{T},
    dotrimming::Bool)::Array{Int64, 1} where {T<:AbstractFloat}

    deleteids, remainingids = _trim(res.integrals)
    @debug "Trimming integration results: $(length(deleteids)) entries out of $(length(res.integrals)) deleted"

    if dotrimming
        deleteat!(res.integrals, deleteids)
        deleteat!(res.volumeID, deleteids)
        res.Y = res.Y[:, remainingids]
    end

    deleteids
end

function split_dataset(dataset::DataSet{T, I})::Tuple{DataSet{T, I}, DataSet{T, I}} where {T<:Real, I<:Integer}
    N = dataset.N
    n = floor(Int64, N / 2)

    ds1 = DataSet(dataset.data[:, 1:n], dataset.logprob[1:n], dataset.weights[1:n], dataset.nsubsets, dataset.subsetsize)
    ds2 = DataSet(dataset.data[:, n+1:N], dataset.logprob[n+1:N], dataset.weights[n+1:N], dataset.nsubsets, dataset.subsetsize)

    return ds1, ds2
end