text:
struct DGModel{BL,G,NFND,NFD,GNF,AS,DS,D,MD}
  balancelaw::BL
  grid::G
  numfluxnondiff::NFND
  numfluxdiff::NFD
  gradnumflux::GNF
  auxstate::AS
  diffstate::DS
  direction::D
  modeldata::MD
end
function DGModel(balancelaw, grid, numfluxnondiff, numfluxdiff, gradnumflux;
                 auxstate=create_auxstate(balancelaw, grid),
                 diffstate=create_diffstate(balancelaw, grid),
                 direction=EveryDirection(), modeldata=nothing)
  DGModel(balancelaw, grid, numfluxnondiff, numfluxdiff, gradnumflux, auxstate,
          diffstate, direction, modeldata)
end

function (dg::DGModel)(dQdt, Q, ::Nothing, t; increment=false)
  bl = dg.balancelaw
  device = typeof(Q.data) <: Array ? CPU() : CUDA()

  grid = dg.grid
  topology = grid.topology

  dim = dimensionality(grid)
  N = polynomialorder(grid)
  Nq = N + 1
  Nqk = dim == 2 ? 1 : Nq
  Nfp = Nq * Nqk
  nrealelem = length(topology.realelems)

  Qvisc = dg.diffstate
  auxstate = dg.auxstate

  FT = eltype(Q)
  nviscstate = num_diffusive(bl, FT)

  lgl_weights_vec = grid.ω
  Dmat = grid.D
  vgeo = grid.vgeo
  sgeo = grid.sgeo
  vmapM = grid.vmapM
  vmapP = grid.vmapP
  elemtobndy = grid.elemtobndy
  polyorder = polynomialorder(dg.grid)

  Np = dofs_per_element(grid)

  communicate = !(isstacked(topology) &&
                  typeof(dg.direction) <: VerticalDirection)

  aux_comm = update_aux!(dg, bl, Q, t)
  @assert typeof(aux_comm) == Bool

  ########################
  # Gradient Computation #
  ########################
  if communicate
    MPIStateArrays.start_ghost_exchange!(Q)
    if aux_comm
      MPIStateArrays.start_ghost_exchange!(auxstate)
    end
  end

  if nviscstate > 0

    @launch(device, threads=(Nq, Nq, Nqk), blocks=nrealelem,
            volumeviscterms!(bl, Val(dim), Val(polyorder), dg.direction, Q.data,
                             Qvisc.data, auxstate.data, vgeo, t, Dmat,
                             topology.realelems))

    if communicate
      MPIStateArrays.finish_ghost_recv!(Q)
      if aux_comm
        MPIStateArrays.finish_ghost_recv!(auxstate)
      end
    end

    @launch(device, threads=Nfp, blocks=nrealelem,
            faceviscterms!(bl, Val(dim), Val(polyorder), dg.direction,
                           dg.gradnumflux, Q.data, Qvisc.data, auxstate.data,
                           vgeo, sgeo, t, vmapM, vmapP, elemtobndy,
                           topology.realelems))

    if communicate
      MPIStateArrays.start_ghost_exchange!(Qvisc)
    end

    aux_comm = update_aux_diffusive!(dg, bl, Q, t)
    @assert typeof(aux_comm) == Bool

    if aux_comm
      MPIStateArrays.start_ghost_exchange!(auxstate)
    end
  end

  ###################
  # RHS Computation #
  ###################
  @launch(device, threads=(Nq, Nq, Nqk), blocks=nrealelem,
          volumerhs!(bl, Val(dim), Val(polyorder), dg.direction, dQdt.data,
                     Q.data, Qvisc.data, auxstate.data, vgeo, t,
                     lgl_weights_vec, Dmat, topology.realelems, increment))

  if communicate
    if nviscstate > 0
      MPIStateArrays.finish_ghost_recv!(Qvisc)
      if aux_comm
        MPIStateArrays.finish_ghost_recv!(auxstate)
      end
    else
      MPIStateArrays.finish_ghost_recv!(Q)
      if aux_comm
        MPIStateArrays.finish_ghost_recv!(auxstate)
      end
    end
  end

  @launch(device, threads=Nfp, blocks=nrealelem,
          facerhs!(bl, Val(dim), Val(polyorder), dg.direction,
                   dg.numfluxnondiff,
                   dg.numfluxdiff,
                   dQdt.data, Q.data, Qvisc.data,
                   auxstate.data, vgeo, sgeo, t, vmapM, vmapP, elemtobndy,
                   topology.realelems))

  # Just to be safe, we wait on the sends we started.
  if communicate
    MPIStateArrays.finish_ghost_send!(Qvisc)
    MPIStateArrays.finish_ghost_send!(Q)
  end
end

function init_ode_state(dg::DGModel, args...;
                        forcecpu=false,
                        commtag=888)
  device = arraytype(dg.grid) <: Array ? CPU() : CUDA()

  bl = dg.balancelaw
  grid = dg.grid

  state = create_state(bl, grid, commtag)

  topology = grid.topology
  Np = dofs_per_element(grid)

  auxstate = dg.auxstate
  dim = dimensionality(grid)
  polyorder = polynomialorder(grid)
  vgeo = grid.vgeo
  nrealelem = length(topology.realelems)

  if !forcecpu
    @launch(device, threads=(Np,), blocks=nrealelem,
            initstate!(bl, Val(dim), Val(polyorder), state.data, auxstate.data, vgeo,
                     topology.realelems, args...))
  else
    h_vgeo = Array(vgeo)
    h_state = similar(state, Array)
    h_auxstate = similar(auxstate, Array)
    h_auxstate .= auxstate
    @launch(CPU(), threads=(Np,), blocks=nrealelem,
      initstate!(bl, Val(dim), Val(polyorder), h_state.data, h_auxstate.data, h_vgeo,
          topology.realelems, args...))
    state .= h_state
  end

  MPIStateArrays.start_ghost_exchange!(state)
  MPIStateArrays.finish_ghost_exchange!(state)

  return state
end

function indefinite_stack_integral!(dg::DGModel, m::BalanceLaw,
                                    Q::MPIStateArray, auxstate::MPIStateArray,
                                    t::Real)

  device = typeof(Q.data) <: Array ? CPU() : CUDA()

  grid = dg.grid
  topology = grid.topology

  dim = dimensionality(grid)
  N = polynomialorder(grid)
  Nq = N + 1
  Nqk = dim == 2 ? 1 : Nq

  FT = eltype(Q)

  vgeo = grid.vgeo
  polyorder = polynomialorder(dg.grid)

  # do integrals
  nintegrals = num_integrals(m, FT)
  nelem = length(topology.elems)
  nvertelem = topology.stacksize
  nhorzelem = div(nelem, nvertelem)

  @launch(device, threads=(Nq, Nqk, 1), blocks=nhorzelem,
          knl_indefinite_stack_integral!(m, Val(dim), Val(polyorder),
                                         Val(nvertelem), Q.data, auxstate.data,
                                         vgeo, grid.Imat, 1:nhorzelem,
                                         Val(nintegrals)))
end

# fallback
function update_aux!(dg::DGModel, bl::BalanceLaw, Q::MPIStateArray, t::Real)
  return false
end

function update_aux_diffusive!(dg::DGModel, bl::BalanceLaw, Q::MPIStateArray, t::Real)
  return false
end

function reverse_indefinite_stack_integral!(dg::DGModel, m::BalanceLaw,
                                            auxstate::MPIStateArray, t::Real)

  device = typeof(auxstate.data) <: Array ? CPU() : CUDA()

  grid = dg.grid
  topology = grid.topology

  dim = dimensionality(grid)
  N = polynomialorder(grid)
  Nq = N + 1
  Nqk = dim == 2 ? 1 : Nq

  FT = eltype(auxstate)

  vgeo = grid.vgeo
  polyorder = polynomialorder(dg.grid)

  # do integrals
  nintegrals = num_integrals(m, FT)
  nelem = length(topology.elems)
  nvertelem = topology.stacksize
  nhorzelem = div(nelem, nvertelem)

  @launch(device, threads=(Nq, Nqk, 1), blocks=nhorzelem,
          knl_reverse_indefinite_stack_integral!(Val(dim), Val(polyorder),
                                                 Val(nvertelem), auxstate.data,
                                                 1:nhorzelem,
                                                 Val(nintegrals)))
end

function nodal_update_aux!(f!, dg::DGModel, m::BalanceLaw, Q::MPIStateArray,
                           t::Real; diffusive=false)
  device = typeof(Q.data) <: Array ? CPU() : CUDA()

  grid = dg.grid
  topology = grid.topology

  dim = dimensionality(grid)
  N = polynomialorder(grid)
  Nq = N + 1
  nrealelem = length(topology.realelems)

  polyorder = polynomialorder(dg.grid)

  Np = dofs_per_element(grid)

  ### update aux variables
  if diffusive
    @launch(device, threads=(Np,), blocks=nrealelem,
            knl_nodal_update_aux!(m, Val(dim), Val(polyorder), f!,
                            Q.data, dg.auxstate.data, dg.diffstate.data, t,
                            topology.realelems))
  else
    @launch(device, threads=(Np,), blocks=nrealelem,
            knl_nodal_update_aux!(m, Val(dim), Val(polyorder), f!,
                            Q.data, dg.auxstate.data, t,
                            topology.realelems))
  end
end

function copy_stack_field_down!(dg::DGModel, m::BalanceLaw,
                                auxstate::MPIStateArray, fldin, fldout)

  device = typeof(auxstate.data) <: Array ? CPU() : CUDA()

  grid = dg.grid
  topology = grid.topology

  dim = dimensionality(grid)
  N = polynomialorder(grid)
  Nq = N + 1
  Nqk = dim == 2 ? 1 : Nq

  DFloat = eltype(auxstate)

  vgeo = grid.vgeo
  polyorder = polynomialorder(dg.grid)

  # do integrals
  nelem = length(topology.elems)
  nvertelem = topology.stacksize
  nhorzelem = div(nelem, nvertelem)

  @launch(device, threads=(Nq, Nqk, 1), blocks=nhorzelem,
          knl_copy_stack_field_down!(Val(dim), Val(polyorder), Val(nvertelem),
                                     auxstate.data, 1:nhorzelem, Val(fldin),
                                     Val(fldout)))
end

function MPIStateArrays.MPIStateArray(dg::DGModel, commtag=888)
  bl = dg.balancelaw
  grid = dg.grid

  state = create_state(bl, grid, commtag)

  return state
end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
export setup_multi_level_PARSDMM

"""
Generate projectors, linear operators, grid information, constraints, set information
on all grids for multilevel PARSDMM
"""
function setup_multi_level_PARSDMM(
                                   m                  ::Vector{TF},
                                   n_levels           ::Integer,
                                   coarsening_factor  ::Union{Real,Integer},
                                   comp_grid,
                                   constraint,
                                   options
                                   ) where {TF<:Real}

TI = Int64

[@spawnat pid comp_grid for pid in workers()] #send computational grid to all workers

#detect 2D or 3D
if length(comp_grid.n)==3 && comp_grid.n[3]>1
  dim3 = true #indicate we are working in 3D
else
  dim3 = false
end

#allocate vectors of vectors and vectors of matrices to store required
#quantities for PARSDMM for each level
m_levels         = Vector{Vector{TF}}(undef,n_levels)
if options.parallel==false
  TD_OP_levels   = Vector{Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}}}(undef,n_levels)
else
  #TD_OP_levels   = Vector{DistributedArrays.DArray{SparseMatrixCSC{TF,TI},1,Array{SparseMatrixCSC{TF,TI},1}}}(n_levels)
  TD_OP_levels   = Vector{DistributedArrays.DArray{Union{JOLI.joAbstractLinearOperator{TF,TF}, SparseMatrixCSC{TF,TI}},1,Array{Union{JOLI.joAbstractLinearOperator{TF,TF}, SparseMatrixCSC{TF,TI}},1}}}(undef,n_levels)
end

AtA_levels       = Vector{Vector{Union{Array{TF, 2},SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}}}(undef,n_levels)
P_sub_levels     = Vector{Vector{Any}}(undef,n_levels)
set_Prop_levels  = Vector{Any}(undef,n_levels)
comp_grid_levels = Vector{Any}(undef,n_levels)

# set up constraints and matrices on coarse level
#1st level is the grid corresponding to the original model m (matrix or tensor)
(P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,TF)
(TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)

# p = length(TD_OP)
# if typeof(AtA)==Vector{Array{TF,2}}
#   AtA_levels       = Vector{Vector{Array{TF,2}}}(undef,n_levels)
# end

i=1
m_levels[i]         = m
TD_OP_levels[i]     = TD_OP
AtA_levels[i]       = AtA
P_sub_levels[i]     = P_sub
set_Prop_levels[i]  = set_Prop
comp_grid_levels[i] = comp_grid

constraint_level = deepcopy(constraint)

for i=2:n_levels

  #new number of grid points on this level
  n = round.( Int, comp_grid.n ./ coarsening_factor^(i-1) )

  #coarsen original model to new grids (Don't need it currently, maybe remove in future version)
  # itp_m     = interpolate(reshape(m,comp_grid.n),BSpline(Linear()))
  #
  # if dim3
  #   m_level = itp_m(range(1,stop=comp_grid.n[1],length=n[1]), range(1,stop=comp_grid.n[2],length=n[2]), range(1,stop=comp_grid.n[3],length=n[3]))
  # else
  #   m_level = itp_m(range(1,stop=comp_grid.n[1],length=n[1]), range(1,stop=comp_grid.n[2],length=n[2]))
  # end
  # m_levels[i] = vec(m_level)

  #create computational grid information for this level
  comp_grid_levels[i]   = deepcopy(comp_grid)
  comp_grid_levels[i].n = n
  comp_grid_levels[i].d = (comp_grid.n ./ n) .* comp_grid.d

  # 'translate' constraints to coarser grids
  # the constraints are originally defined on the grid of the model on the finest grid.
  # We need to define the equivalent constraints on coarser grids
  constraint_level = constraint2coarse(constraint_level,comp_grid_levels[i],coarsening_factor)

  #set up constraints on new level
  println(TF)
  (P_sub_l,TD_OP_l,set_Prop_l)  = setup_constraints(constraint_level,comp_grid_levels[i],TF)
  (TD_OP_l,AtA_l,dummy1,dummy2) = PARSDMM_precompute_distribute(TD_OP_l,set_Prop_l,comp_grid_levels[i],options)

  #save information
  TD_OP_levels[i]    = TD_OP_l
  AtA_levels[i]      = AtA_l
  P_sub_levels[i]    = P_sub_l
  set_Prop_levels[i] = set_Prop_l

end #for loop

# #replicate edges of the subsampled m to avoid artifacts from the blur kernel
# for i=1:n_levels
#   m_levels[i]=reshape(m_levels[i],comp_grid_levels[i].n)
#   if dim3 == true
#     m_levels[i][1,:,:]=m_levels[i][2,:,:]
#     m_levels[i][1,:,:]=m_levels[i][2,:,:]
#     m_levels[i][end,:,:]=m_levels[i][end-1,:,:]
#     m_levels[i][end,:,:]=m_levels[i][end-1,:,:]
#
#     m_levels[i][:,1,:]=m_levels[i][:,2,:]
#     m_levels[i][:,1,:]=m_levels[i][:,2,:]
#     m_levels[i][:,end,:]=m_levels[i][:,end-1,:]
#     m_levels[i][:,end,:]=m_levels[i][:,end-1,:]
#
#     m_levels[i][:,:,1]=m_levels[i][:,:,2]
#     m_levels[i][:,:,1]=m_levels[i][:,:,2]
#     m_levels[i][:,:,end]=m_levels[i][:,:,end-1]
#     m_levels[i][:,:,end]=m_levels[i][:,:,end-1]
#   else
#     m_levels[i][1,:]=m_levels[i][2,:]
#     m_levels[i][1,:]=m_levels[i][2,:]
#     m_levels[i][end,:]=m_levels[i][end-1,:]
#     m_levels[i][end,:]=m_levels[i][end-1,:]
#
#     m_levels[i][:,1]=m_levels[i][:,2]
#     m_levels[i][:,1]=m_levels[i][:,2]
#     m_levels[i][:,end]=m_levels[i][:,end-1]
#     m_levels[i][:,end]=m_levels[i][:,end-1]
#   end
#   m_levels[i]=vec(m_levels[i])
# end

#return m_levels, TD_OP_levels, AtA_levels, P_sub_levels, set_Prop_levels, comp_grid_levels, constraint_level
return TD_OP_levels, AtA_levels, P_sub_levels, set_Prop_levels, comp_grid_levels, constraint_level

end #end setup_multi_level_PARSDMM

ext: jl

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
using StatsBase, Distributions, Random, LinearAlgebra
Random.seed!(1)

function crudeSimulation(deltaT,T,Q,initProb)
    n = size(Q)[1]
    Pdelta = I + Q*deltaT
    state  = sample(1:n,weights(initProb))
    t = 0.0
    while t < T
        t += deltaT
        state = sample(1:n,weights(Pdelta[state,:]))
    end
    return state
end

function doobGillespie(T,Q,initProb)
    n = size(Q)[1]
    Pjump  = (Q-diagm(0 => diag(Q)))./-diag(Q)
    lamVec = -diag(Q)
    state  = sample(1:n,weights(initProb))
    sojournTime = rand(Exponential(1/lamVec[state]))
    t = 0.0
    while t + sojournTime < T
        t += sojournTime
        state = sample(1:n,weights(Pjump[state,:]))
        sojournTime = rand(Exponential(1/lamVec[state]))
    end
    return state
end

T, N = 0.25, 10^5

Q = [-3 1 2
     1 -2 1
     0 1.5 -1.5]

p0 = [0.4 0.5 0.1]

crudeSimEst = counts([crudeSimulation(10^-3., T, Q, p0) for _ in 1:N])/N
doobGillespieEst = counts([doobGillespie(T, Q, p0) for _ in 1:N])/N
explicitEst = p0*exp(Q*T)

println("CrudeSim: \t\t", crudeSimEst)
println("Doob Gillespie Sim: \t", doobGillespieEst)
println("Explicit: \t\t", explicitEst)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.DGmethods
using CLIMA.DGmethods.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using CLIMA.Atmos
using CLIMA.VariableTemplates
using CLIMA.MoistThermodynamics
using CLIMA.PlanetParameters
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates
using CLIMA.VTK

if !@isdefined integration_testing
  const integration_testing =
    parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
end

using CLIMA.Atmos
using CLIMA.Atmos: internal_energy, thermo_state
import CLIMA.Atmos: MoistureModel, temperature, pressure, soundspeed

init_state!(bl, state, aux, coords, t) = nothing

# initial condition
using CLIMA.Atmos: vars_aux

function run1(mpicomm, ArrayType, dim, topl, N, timeend, FT, dt)

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N
                                         )

  T_s = 320.0
  RH = 0.01
  model = AtmosModel{FT}(AtmosLESConfiguration;
                         ref_state=HydrostaticState(IsothermalProfile(T_s), RH),
                        init_state=init_state!)

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

  Q = init_ode_state(dg, FT(0))

  mkpath("vtk")
  outprefix = @sprintf("vtk/refstate")
  writevtk(outprefix, dg.auxstate, dg, flattenednames(vars_aux(model, FT)))
  return FT(0)
end

function run2(mpicomm, ArrayType, dim, topl, N, timeend, FT, dt)

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N
                                         )

  T_min, T_s, Γ = FT(290), FT(320), FT(6.5*10^-3)
  RH = 0.01
  model = AtmosModel{FT}(AtmosLESConfiguration;
                         ref_state=HydrostaticState(LinearTemperatureProfile(T_min, T_s, Γ), RH),
                         init_state=init_state!)

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

  Q = init_ode_state(dg, FT(0))

  mkpath("vtk")
  outprefix = @sprintf("vtk/refstate")
  writevtk(outprefix, dg.auxstate, dg, flattenednames(vars_aux(model, FT)))
  return FT(0)
end

using Test
let
  CLIMA.init()
  ArrayType = CLIMA.array_type()

  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
    ll == "WARN"  ? Logging.Warn  :
    ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))

  polynomialorder = 4
  base_num_elem = 4

  expected_result = [0.0 0.0 0.0; 0.0 0.0 0.0]
  lvls = integration_testing ? size(expected_result, 2) : 1

for FT in (Float64,) #Float32)
  result = zeros(FT, lvls)
  x_max = FT(25*10^3)
  y_max = FT(25*10^3)
  z_max = FT(25*10^3)
  dim = 3
  for l = 1:lvls
    Ne = (2^(l-1) * base_num_elem, 2^(l-1) * base_num_elem)
    brickrange = (range(FT(0); length=Ne[1]+1, stop=x_max),
                  range(FT(0); length=Ne[2]+1, stop=y_max),
                  range(FT(0); length=Ne[2]+1, stop=z_max))
    topl = BrickTopology(mpicomm, brickrange,
                         periodicity = (false, false, false))
    dt = 5e-3 / Ne[1]

    timeend = 2*dt
    nsteps = ceil(Int64, timeend / dt)
    dt = timeend / nsteps

    @info (ArrayType, FT, dim)
    result[l] = run1(mpicomm, ArrayType, dim, topl,
                    polynomialorder, timeend, FT, dt)
    result[l] = run2(mpicomm, ArrayType, dim, topl,
                    polynomialorder, timeend, FT, dt)
  end
  if integration_testing
    @info begin
      msg = ""
      for l = 1:lvls-1
        rate = log2(result[l]) - log2(result[l+1])
        msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
      end
      msg
    end
  end
end
end

#nothing

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
using Test, SafeTestsets

@time @safetestset "Unit Tests" begin include("unit_tests/unit_tests.jl") end
@time @safetestset "Example Models" begin include("examples/examples.jl") end

if ENV["JULIA_EXTENDED_TESTS"] == "true"

end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using GaloisFields, LinearAlgebra
const F5 = @GaloisField 5
A = F5[1 2; 0 3]

# %%
det(A)

# %%
B = F5[1 2; 3 4]

# %%
det(B)

# %%
F5(3) |> dump

# %%
F5 |> dump

# %%
Base.abs(x::GaloisFields.PrimeField) = x.n # type-piracy!
det(B)

# %%
F3 = @GaloisField 3
[det(F3.(rand(0:2, 3, 3))) for _ in 1:20]

# %%
Base.abs(x::GaloisFields.ExtensionField) = x.coeffs # type-piracy!
Base.isless(x::GaloisFields.PrimeField, y::GaloisFields.PrimeField) = isless(x.n, y.n) # type-piracy!
F9 = @GaloisField! 3 a^2 + 1
C = F9[1+a 2+a; 2+a 2+a]

# %%
det(C)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using ConstrainedDynamics

# Body
Box(rand(4)...; color = RGBA(rand(3)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
@test true
Cylinder(rand(3)...; color = RGBA(rand(3)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
@test true
Sphere(rand(2)...; color = RGBA(rand(3)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
@test true
Pyramid(rand(3)...; color = RGBA(rand(3)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
@test true
Mesh("test.obj", rand(), rand(3,3); color = RGBA(rand(3)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
@test true
# Just Shape
Box(rand(3)...; color = RGBA(rand(3)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
@test true
Cylinder(rand(2)...; color = RGBA(rand(3)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
@test true
Sphere(rand()...; color = RGBA(rand(3)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
@test true
Pyramid(rand(2)...; color = RGBA(rand(3)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
@test true
Mesh("test.obj"; color = RGBA(rand(3)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
@test true

storage1 = Storage{Float64}(Base.OneTo(2),2)
@test true
storage2 = Storage([[[[0.0;0;0] for i=1:2]];[[[0.0;0;0] for i=1:2]]],[[[UnitQuaternion(RotX(0.0)) for i=1:2]];[[UnitQuaternion(RotX(0.0)) for i=1:2]]])
@test true


text:
# * Linear interpolation

# https://github.com/JuliaLang/julia/pull/18777
lerp(a::T,b::T,t) where T = T(fma(t, b, fma(-t, a, a)))
lerp(a::R,b::R,t::C) where {R<:Real,C<:Complex} = lerp(a,b,real(t)) + im*lerp(a,b,imag(t))
lerp(a::C,b::C,t::R) where {R<:Real,C<:Complex} = lerp(real(a),real(b),t) + im*lerp(imag(a),imag(b),(t))

# * Gaußian quadrature

"""
    change_interval!(xs, ws, x, w[, a=0, b=1, γ=1])

Transform the Gaußian quadrature roots `x` and weights `w` on the
elementary interval `[-1,1]` to the interval `[γ*a,γ*b]` and store the
result in `xs` and `ws`, respectively. `γ` is an optional root of
unity, used to complex-rotate the roots (but not the weights).
"""
function change_interval!(xs::AbstractVector{T}, ws, x, w,
                          a=zero(T), b=one(T), γ=one(T)) where T
    xs .= lerp.(γ*a, γ*b, (x .+ 1)/2)
    ws .= (b-a)*w/2
    xs,ws
end

change_interval(x::AbstractVector{T}, w::AbstractVector{T},
                a=zero(T), b=one(T), γ::U=one(T)) where {T,U} =
    change_interval!(similar(x,U), similar(w), x, w, a, b)

# * Gauß–Legendre quadrature

"""
    num_quadrature_points(k, k′)

The number of quadrature points needed to exactly compute the matrix
elements of an operator of polynomial order `k′` with respect to a
basis of order `k`.
"""
function num_quadrature_points(k, k′)
    N2 = 2*(k-1) + k′
    N2>>1 + N2&1
end

"""
    lgwt(t, N) -> (x,w)

Generate the `N` Gauß–Legendre quadrature roots `x` and associated
weights `w`, with respect to the B-spline basis generated by the knot
set `t`.

# Examples

```jldoctest
julia> CompactBases.lgwt(LinearKnotSet(2, 0, 1, 3), 2)
([0.0704416, 0.262892, 0.403775, 0.596225, 0.737108, 0.929558], [0.166667, 0.166667, 0.166667, 0.166667, 0.166667, 0.166667])

julia> CompactBases.lgwt(ExpKnotSet(2, -4, 2, 7), 2)
([2.11325e-5, 7.88675e-5, 0.000290192, 0.000809808, 0.00290192, 0.00809808, 0.0290192, 0.0809808, 0.290192, 0.809808, 2.90192, 8.09808, 29.0192, 80.9808], [5.0e-5, 5.0e-5, 0.00045, 0.00045, 0.0045, 0.0045, 0.045, 0.045, 0.45, 0.45, 4.5, 4.5, 45.0, 45.0])
```
"""
function lgwt(t::AbstractKnotSet{k,ml,mr,T}, N) where {k,ml,mr,T}
    2N-1 ≥ 2(k-1) || @warn "N = $N quadrature point$(N > 1 ? "s" : "") not enough to calculate overlaps between polynomials of order k = $k"
    x, w = gausslegendre(N)

    nei = nonempty_intervals(t)
    ni = length(nei)
    xo = zeros(T, ni*length(x))
    wo = zeros(T, ni*length(x))

    for (i,j) in enumerate(nei)
        sel = (i-1)*N+1 : i*N
        change_interval!(view(xo, sel), view(wo, sel),
                         x, w,
                         t[j], t[j+1])
    end

    xo,wo
end

# * Gauß–Lobatto

function element_grid(order, a::T, b::T, c::T=zero(T), eiϕ=one(T)) where T
    x,w = gausslobatto(order)
    xs,ws = change_interval(x, w, a-c, b-c, eiϕ)
    c .+ xs, ws
end


text:
module CodecZstd

export
    ZstdCompressor,
    ZstdCompressorStream,
    ZstdDecompressor,
    ZstdDecompressorStream

import TranscodingStreams:
    TranscodingStreams,
    TranscodingStream,
    Memory,
    Error,
    initialize,
    finalize,
    splitkwargs

using Zstd_jll

include(joinpath("libzstd","LibZstd.jl"))
include("libzstd.jl")
include("compression.jl")
include("decompression.jl")

end # module


text:
module BoundaryTriangulationsTests

using Test
using Gridap.Helpers
using Gridap.Fields
using Gridap.Arrays
using Gridap.ReferenceFEs
using Gridap.Geometry

domain = (0,4,0,4)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

btrian = BoundaryTriangulation(model)
test_boundary_triangulation(btrian)

cellids = collect(1:num_cells(model))

face_to_cellid = reindex(cellids,btrian)
@test face_to_cellid == get_face_to_cell(btrian)

trian = get_volume_triangulation(btrian)
q2x = get_cell_map(trian)
s2x = restrict(q2x,btrian)

s = CompressedArray([Point{1,Float64}[(0.25,),(0.75,)]],get_cell_type(btrian))
x = evaluate(s2x,s)
r = Vector{Point{2,Float64}}[
  [(0.5,0.0),(1.5,0.0)],[(0.0,0.5),(0.0,1.5)],
  [(2.5,0.0),(3.5,0.0)],[(4.0,0.5),(4.0,1.5)],
  [(0.5,4.0),(1.5,4.0)],[(0.0,2.5),(0.0,3.5)],
  [(2.5,4.0),(3.5,4.0)],[(4.0,2.5),(4.0,3.5)]]
test_array(x,r)

nvec = get_normal_vector(btrian)
nvec_x = evaluate(nvec,s)

domain = (0,4,0,4)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

btrian = BoundaryTriangulation(model,"tag_8")
test_boundary_triangulation(btrian)
@test btrian.face_trian.cell_to_oldcell == get_face_to_face(btrian)
@test 1 == get_cell_around(btrian)

btrian = BoundaryTriangulation(model,get_face_labeling(model),"tag_8")
test_boundary_triangulation(btrian)

s = CompressedArray([Point{1,Float64}[(0.25,),(0.75,)]],get_cell_type(btrian))
nvec = get_normal_vector(btrian)
nvec_x = evaluate(nvec,s)
s2x = get_cell_map(btrian)
x = evaluate(s2x,s)

@test get_cell_id(btrian) == get_face_to_cell(btrian)
r = rand(num_cells(trian))
@test reindex(r,btrian) == r[get_face_to_cell(btrian)]

end # module


text:
using YaoLang
using Test
using FFTW
using YaoArrayRegister

@testset "runtime" begin
    include("runtime/locations.jl")
end

@testset "compiler" begin
    include("compiler/circuit.jl")
    include("compiler/validation.jl")
end


text:
const NGROUP       = 6
const MAXOVERLAY   = 500

const MAXLINE      = 100
const MAXLINEPNT   = 1000
const MAXPLANEGRID = 200

@enum mjtCatBit begin           # bitflags for vGeom category
    CAT_STATIC        = 1  # model elements in body 0
    CAT_DYNAMIC       = 2  # model elements in all other bodies
    CAT_DECOR         = 4  # decorative geoms
    CAT_ALL           = 7  # select all categories
end

@enum mjtMouse begin            # mouse interaction mode
    MOUSE_NONE         = 0 # no action
    MOUSE_ROTATE_V         # rotate, vertical plane
    MOUSE_ROTATE_H         # rotate, horizontal plane
    MOUSE_MOVE_V           # move, vertical plane
    MOUSE_MOVE_H           # move, horizontal plane
    MOUSE_ZOOM             # zoom
    MOUSE_SELECT           # selection
end

@enum mjtPertBit begin          # mouse perturbations
    PERT_TRANSLATE    = 1  # translation
    PERT_ROTATE       = 2  # rotation
end

@enum mjtCamera begin           # abstract camera type
    CAMERA_FREE        = 0 # free camera
    CAMERA_TRACKING        # tracking camera; uses trackbodyid
    CAMERA_FIXED           # fixed camera; uses fixedcamid
    CAMERA_USER            # user is responsible for setting OpenGL camera
end

@enum mjtLabel begin            # object labeling
    LABEL_NONE        = 0  # nothing
    LABEL_BODY             # body labels
    LABEL_JOINT            # joint labels
    LABEL_GEOM             # geom labels
    LABEL_SITE             # site labels
    LABEL_CAMERA           # camera labels
    LABEL_LIGHT            # light labels
    LABEL_TENDON           # tendon labels
    LABEL_ACTUATOR         # actuator labels
    LABEL_CONSTRAINT       # constraint labels
    LABEL_SKIN             # skin labels
    LABEL_SELECTION        # selected object
    LABEL_SELPNT           # coordinates of selection point
    LABEL_CONTACTFORCE     # magnitude of contact force

    NLABEL                 # number of label types
end

@enum mjtFrame begin            # frame visualization
    FRAME_NONE        = 0  # no frames
    FRAME_BODY             # body frames
    FRAME_GEOM             # geom frames
    FRAME_SITE             # site frames
    FRAME_CAMERA           # camera frames
    FRAME_LIGHT            # light frames
    FRAME_WORLD            # world frame

    NFRAME                 # number of visualization frames
end

@enum mjtVisFlag begin          # flags enabling model element visualization
    VIS_CONVEXHULL    = 0  # mesh convex hull
    VIS_TEXTURE            # textures
    VIS_JOINT              # joints
    VIS_ACTUATOR           # actuators
    VIS_CAMERA             # cameras
    VIS_LIGHT              # lights
    VIS_TENDON             # tendons
    VIS_RANGEFINDER        # rangefinder sensors
    VIS_CONSTRAINT         # point constraints
    VIS_INERTIA            # equivalent inertia boxes
    VIS_SCLINERTIA         # scale equivalent inertia boxes with mass
    VIS_PERTFORCE          # perturbation force
    VIS_PERTOBJ            # perturbation object
    VIS_CONTACTPOINT       # contact points
    VIS_CONTACTFORCE       # contact force
    VIS_CONTACTSPLIT       # split contact force into normal and tanget
    VIS_TRANSPARENT        # make dynamic geoms more transparent
    VIS_AUTOCONNECT        # auto connect joints and body coms
    VIS_COM                # center of mass
    VIS_SELECT             # selection point
    VIS_STATIC             # static bodies
    VIS_SKIN               # skin

    NVISFLAG               # number of visualization flags
end

@enum mjtRndFlag begin          # flags enabling rendering effects
    RND_SHADOW        = 0  # shadows
    RND_WIREFRAME          # wireframe
    RND_REFLECTION         # reflections
    RND_ADDITIVE           # additive transparency
    RND_SKYBOX             # skybox
    RND_FOG                # fog
    RND_HAZE               # haze
    RND_SEGMENT            # segmentation with random color
    RND_IDCOLOR            # segmentation with segid color

    NRNDFLAG               # number of rendering flags
end

@enum mjtStereo begin           # type of stereo rendering
    STEREO_NONE       = 0  # no stereo; use left eye only
    STEREO_QUADBUFFERED    # quad buffered; revert to side-by-side if no hardware support
    STEREO_SIDEBYSIDE      # side-by-side
end

struct mjvPerturb
   select::Cint
   skinselect::Cint
   active::Cint
   refpos::SVector{3, mjtNum}
   refquat::SVector{4, mjtNum}
   localpos::SVector{3, mjtNum}
   scale::mjtNum
end

struct mjvCamera
   _type::Cint
   fixedcamid::Cint
   trackbodyid::Cint
   lookat::SVector{3, mjtNum}
   distance::mjtNum
   azimuth::mjtNum
   elevation::mjtNum
end

struct mjvGLCamera
   pos::SVector{3, Cfloat}
   forward::SVector{3, Cfloat}
   up::SVector{3, Cfloat}
   frustum_center::Cfloat
   frustum_bottom::Cfloat
   frustum_top::Cfloat
   frustum_near::Cfloat
   frustum_far::Cfloat
end

struct mjvGeom
   _type::Cint
   dataid::Cint
   objtype::Cint
   objid::Cint
   category::Cint
   texid::Cint
   texuniform::Cint
   texcoord::Cint
   segid::Cint
   texrepeat::SVector{2, Cfloat}
   size::SVector{3, Cfloat}
   pos::SVector{3, Cfloat}
   mat::SVector{9, Cfloat}
   rgba::SVector{4, Cfloat}
   emission::Cfloat
   specular::Cfloat
   shininess::Cfloat
   reflectance::Cfloat
   label::SVector{100, UInt8}
   camdist::Cfloat
   rbound::Cfloat
   transparent::mjtByte
end

struct mjvLight
   pos::SVector{3, Cfloat}
   dir::SVector{3, Cfloat}
   attenuation::SVector{3, Cfloat}
   cutoff::Cfloat
   exponent::Cfloat
   ambient::SVector{3, Cfloat}
   diffuse::SVector{3, Cfloat}
   specular::SVector{3, Cfloat}
   headlight::mjtByte
   directional::mjtByte
   castshadow::mjtByte
end

struct mjvOption
   label::Cint
   frame::Cint
   geomgroup::SVector{NGROUP, mjtByte}
   sitegroup::SVector{NGROUP, mjtByte}
   jointgroup::SVector{NGROUP, mjtByte}
   tendongroup::SVector{NGROUP, mjtByte}
   actuatorgroup::SVector{NGROUP, mjtByte}
   flags::SVector{Int(NVISFLAG), mjtByte}
end

struct mjvScene
   maxgeom::Cint
   ngeom::Cint
   geoms::Ptr{mjvGeom}
   geomorder::Ptr{Cint}

   nskin::Cint
   skinfacenum::Ptr{Cint}
   skinvertadr::Ptr{Cint}
   skinvertnum::Ptr{Cint}
   skinvert::Ptr{Cfloat}
   skinnormal::Ptr{Cfloat}

   nlight::Cint
   lights::SVector{8, mjvLight}
   camera::SVector{2, mjvGLCamera}
   enabletransform::mjtByte
   translate::SVector{3, Cfloat}
   rotate::SVector{4, Cfloat}
   scale::Cfloat
   stereo::Cint
   flags::SVector{5, mjtByte}
end

struct mjvFigure
   flg_legend::Cint
   flg_ticklabel::SVector{2, Cint}
   flg_extend::Cint
   flg_barplot::Cint
   flg_selection::Cint
   flg_symmetric::Cint

   legendoff::Cint
   gridsize::SVector{2, Cint}
   selection::Cint
   highlight::SVector{2, Cint}
   gridrgb::SVector{3, Cfloat}
   gridwidth::Cfloat
   figurergba::SVector{4, Cfloat}
   panergba::SVector{4, Cfloat}
   legendrgba::SVector{4, Cfloat}
   textrgb::SVector{3, Cfloat}
	range::SVector{2, SVector{2, Cfloat}}
   xlabel::SVector{100, UInt8}
   title::SVector{100, UInt8}
   xformat::SVector{20, UInt8}
   yformat::SVector{20, UInt8}
   minwidth::SVector{20, UInt8}

   linepnt::SVector{MAXLINE, Cint}
	linergb::SVector{MAXLINE, SVector{3, Cfloat}}
   linewidth::SVector{MAXLINE, Cfloat}
	linedata::SVector{MAXLINE, SVector{2*MAXLINEPNT, Cfloat}}
	linename::SVector{MAXLINE, SVector{100, UInt8}}

   xaxispixel::SVector{2, Cint}
   yaxispixel::SVector{2, Cint}
   xaxisdata::SVector{2, Cint}
   yaxisdata::SVector{2, Cint}
end


text:
sig = [cos(2*pi*0.15*i) for i=1:2048];
@test_approx_eq_eps Estimation.quinn_fernandes(sig,2*pi*0.1)/(2*pi) 0.15 1e-5


text:
mutable struct BatchEncoder <: SEALObject
  handle::Ptr{Cvoid}
  context::SEALContext

  function BatchEncoder(context)
    handleref = Ref{Ptr{Cvoid}}(C_NULL)
    retval = ccall((:BatchEncoder_Create, libsealc), Clong,
                   (Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                   context, handleref)
    @check_return_value retval
    return BatchEncoder(handleref[], context)
  end

  function BatchEncoder(handle::Ptr{Cvoid}, context)
    object = new(handle, context)
    finalizer(destroy!, object)
    return object
  end
end

function destroy!(object::BatchEncoder)
  if isallocated(object)
    @check_return_value ccall((:BatchEncoder_Destroy, libsealc), Clong, (Ptr{Cvoid},), object)
    sethandle!(object, C_NULL)
  end

  return nothing
end

function slot_count(encoder::BatchEncoder)
  count = Ref{UInt64}(0)
  retval = ccall((:BatchEncoder_GetSlotCount, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}),
                 encoder, count)
  @check_return_value retval
  return Int(count[])
end

function encode!(destination::Plaintext, values::DenseArray{UInt64}, encoder::BatchEncoder)
  retval = ccall((:BatchEncoder_Encode1, libsealc), Clong,
                 (Ptr{Cvoid}, UInt64, Ref{UInt64}, Ptr{Cvoid}),
                 encoder, length(values), values, destination)
  @check_return_value retval
  return destination
end

function encode!(destination::Plaintext, values::DenseArray{Int64}, encoder::BatchEncoder)
  retval = ccall((:BatchEncoder_Encode2, libsealc), Clong,
                 (Ptr{Cvoid}, UInt64, Ref{Int64}, Ptr{Cvoid}),
                 encoder, length(values), values, destination)
  @check_return_value retval
  return destination
end

function decode!(destination::DenseVector{UInt64}, plain::Plaintext, encoder::BatchEncoder)
  count = Ref{UInt64}(0)
  retval = ccall((:BatchEncoder_Decode1, libsealc), Clong,
                 (Ptr{Cvoid}, Ptr{Cvoid}, Ref{UInt64}, Ref{UInt64}, Ptr{Cvoid}),
                 encoder, plain, count, destination, C_NULL)
  @check_return_value retval
  return destination
end

function decode!(destination::DenseVector{Int64}, plain::Plaintext, encoder::BatchEncoder)
  count = Ref{UInt64}(0)
  retval = ccall((:BatchEncoder_Decode2, libsealc), Clong,
                 (Ptr{Cvoid}, Ptr{Cvoid}, Ref{UInt64}, Ref{Int64}, Ptr{Cvoid}),
                 encoder, plain, count, destination, C_NULL)
  @check_return_value retval
  return destination
end


text:
using HorizonSideRobots
include("Functions.jl")

function paint_stairs!(r::Robot)
    o = move_with_path!(r, Ost)
    s = move_with_path!(r, Sud)
    len = move_with_path!(r, West)
    len_cur_step = len
    num_side = 3

    while (!isborder(r, Nord)) && (len_cur_step > 0)
        move_n_steps_putmarkers!(r, HorizonSide(num_side), len_cur_step)
        len_cur_step -= 1
        move!(r, Nord)
        num_side = (num_side + 2) % 4
        if num_side == 1
            move!(r, West)
        end
    end

    if isborder(r, Nord)
        if num_side == 1
            move!(r, West)
        end
        move_n_steps_putmarkers!(r, HorizonSide(num_side), len_cur_step)
    end

    movements!(r, West)
    movements!(r, Sud)
    movements!(r, Ost, len - o)
    movements!(r, Nord, s)
end

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
type solution
    seq::Array{Int64}
    fitness::Float64
end

function TwoPointCrossOver(best::Array{Int64,1},sol::Array{Int64,1})
        p_1=rand(2:length(best)-1)
        p_2=rand(2:length(best)-1)

        if p_1==p_2
            p_2=p_1+1
        end
        int_best=best[p_1:p_2]

        new_sol=deepcopy(sol)
        new_sol[p_1:p_2]=int_best
        return new_sol
end

function Sort(population) # sort population based of their fitness
    fits=zeros(length(population))
    for i=1:length(population)
        fits[i]=population[i].fitness
    end
    key=sortperm(fits)
    s_population=population[key];
end

function Partition(s_population,m,n) # partition population into m memplex
    p_population=[]
    memplex=[];
    s=1
    c=1
    for j=1:m
        s=c
        while s<=m*n
            push!(memplex,s_population[s])
            s=s+m
        end
        push!(p_population,memplex)
        memplex=[];
        c=c+1
    end
    return p_population
end

function HammingDistance(a::Array{Int32},b::Array{Int32})
	distance=countnz(a-b);
end

function FrogLeap(best_sol::Array{Int64},worst_sol::Array{Int64},smax::Int64)
	new_sol=deepcopy(worst_sol)
	p1=rand(2:length(best_sol)-smax)
	p2=p1+min(HammingDistance(best_sol,worst_sol),smax)

	@show block=best_sol[p1:p2];
	@show new_sol[p1:p2]=block
	return new_sol
end

function SelectionProb(n)
    prob_vector=[];
    for i=1:n
        p=2*(n+1-i)/(n*(n+1))
        push!(prob_vector,p)
    end
    return prob_vector
end

function CreatePop(m,n,K)
    population=[];
    for i=1:n*m
        seq=zeros(K);
        for j=1:K
            if rand()<0.5
               seq[j]=1
           else
               seq[j]=0
            end
        end
        sol=solution(seq,0)
        push!(population,sol)
    end
    return population
end

max_stars_repo_name: imanpanahi/MilProject

max_issues_repo_name: imanpanahi/MilProject

max_forks_repo_name: imanpanahi/MilProject


text:
function JacobianPropagatorConverter(arg0::NumericalPropagatorBuilder, arg1::jdouble, arg2::jint)
    return JacobianPropagatorConverter((NumericalPropagatorBuilder, jdouble, jint), arg0, arg1, arg2)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule GAP_pkg_orb_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("GAP_pkg_orb")
JLLWrappers.@generate_main_file("GAP_pkg_orb", UUID("b689cd95-dcfd-544e-a174-ce2fd787c795"))
end  # module GAP_pkg_orb_jll


text:
#script for making visualizations

function visualize_observations(objects_observed::Matrix{Array{Array{Detection2D}}}, v::Int64, f::Int64, receptive_fields)

    xs = []
    ys = []
    cs = []
    for rf = 1:length(receptive_fields)
        objects = objects_observed[v, f][rf]
        if length(objects) > 0
            for i = 1:length(objects)
                obj = objects[i]
                push!(xs, obj[1])
                push!(ys, obj[2])
                push!(cs, obj[3])
            end
        end
    end

    if length(xs) < 1
        push!(xs, Inf)
        push!(ys, Inf)
        push!(cs, 1000000)
    end

    println("xs", xs)
    println("ys", ys)
    pyplot()
    p = scatter(xs, ys, color = cs, series_annotations = text.(cs, :bottom), title = "Observations/Detections", xlims = (0, 360), ylims = (0, 240))
    display(p)
    savefig(p, "observations.pdf")
end

#j is which trace
function visualize_trace(traces, j::Int64, camera_trajectories::Matrix{Camera_Params}, v::Int64, f::Int64, params::Video_Params)
    inferred_scene = traces[j][:videos => v => :init_scene]
    println("inferred_scene ", inferred_scene)

    #figure out where each object in the inferred_scene scene would show up in frame f
    xs = []
    ys = []
    cs = []
    if length(inferred_scene) > 0
        for i = 1:length(inferred_scene)
            obj = inferred_scene[i]
            (x,y) = get_image_xy(camera_trajectories[v, f], params, Coordinate(obj[1], obj[2], obj[3]))
            push!(xs, x)
            push!(ys, y)
            push!(cs, obj[4])
        end
    else
        push!(xs, Inf)
        push!(ys, Inf)
        push!(cs, 1000000)
    end

    #how many objects are out-of-view? total - number objects in-view
    in_x = (xs .<= 360) .& (xs .>= 0)
    in_y = (ys .<= 240) .& (ys .>= 0)
    num_obj_out_of_view = length(inferred_scene) - sum(in_x .& in_y)

    p = plt.figure()
    plt.scatter(xs, ys)
    for (i, c) in enumerate(cs)
        plt.annotate(c, (xs[i], ys[i]))
    end
    plt.xlim([0, 360])
    plt.ylim([0, 240])
    p.suptitle("particle $(j)'s inferences. weight: $(get_score(traces[j])). num obj out of view: $num_obj_out_of_view")

    #add observations. hard-coded.
    # circle1 = plt.Circle((160, 120), 30, color='r')
    # circle2 = plt.Circle((165, 125), 30, color='b')
    # fig = plt.gcf()
    # ax = fig.gca()
    # ax.add_patch(circle1)
    # ax.add_patch(circle2)
    # println("here")

    p.savefig("v $(v) particle $(j) .pdf")
end

function visualize_trace_with_heatmap(traces, camera_trajectories::Matrix{Camera_Params}, v::Int64, f::Int64, params::Video_Params)
    xs = []
    ys = []
    cs = []
    weights =[]

    for (j, trace) in enumerate(traces)
        inferred_scene = trace[:videos => v => :init_scene]
        for obj in inferred_scene
            (x, y) = get_image_xy(camera_trajectories[v, f], params, Coordinate(obj[1], obj[2], obj[3]))
            push!(xs, x)
            push!(ys, y)
            push!(cs, obj[4])
            push!(weights, get_score(trace))
        end
        visualize_scene_3d(inferred_scene, "inference_particle_$(j)", "inference_particle_$(j).pdf")
    end

    # Makes heat map that ignores object categories
    p = plt.figure()
    plt.xlim([0, 360])
    plt.ylim([0, 240])
    bins = plt.hexbin(xs, ys, C=weights, gridsize=30, cmap="Greens")
    plt.suptitle("Video $(v)")
    p.savefig("heatmap.pdf")
    plt.close()

    # Makes heat map for each object category
    for cat in Set(cs)
        cat_xs = [x for (x, c) in zip(xs, cs) if c==cat]
        cat_ys = [y for (y, c) in zip(ys, cs) if c==cat]
        cat_weights = [weight for (weight, c) in zip(weights, cs) if c==cat]

        q = plt.figure()
        plt.xlim([0, 360])
        plt.ylim([0, 240])
        plt.hexbin(cat_xs, cat_ys, C=cat_weights, gridsize=10, cmap="Greens")
        q.suptitle("Video $(v) Object $(cat)")
        q.savefig("heatmap_cat$(cat).pdf")
        plt.close()
    end
end

function visualize_scene_3d(scene, title, filename)
    xs = []
    ys = []
    zs = []
    cs = []

    for obj in scene
        push!(xs, obj[1])
        push!(ys, obj[2])
        push!(zs, obj[3])
        push!(cs, obj[4])
    end

    p = plt.figure()
    ax = p.add_subplot(projection="3d")
    ax.scatter(xs, ys, zs)
    ax.set_xlim(-16, 10)
    ax.set_ylim(-10, 10)
    ax.set_zlim(-10, 10)

    for (i, c) in enumerate(cs)
        ax.text(xs[i], ys[i], zs[i], c)
    end
    ax.set_title(title)

    p.savefig(filename)
end

export visualize_observations
export visualize_trace
export visualize_trace_with_heatmap
export visualize_scene_3d

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
#=――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
    NOTE
    The algorithm used for inflecting English words into plurals is descibed by
    Damian Conway of Monash University. The original paper is accessible at
    <http://www.csse.monash.edu.au/~damian/papers/HTML/Plurals.html>

    Several of the tables desribed within are obtained from
    <http://www.csse.monash.edu.au/~damian/papers/HTML/Plurals_AppendixA.html>

    CITATION
    Conway, D. M. (1998, August). An algorithmic approach to english
    pluralization. In Proceedings of the Second Annual Perl Conference.
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――=#

const SUFFIX_NOINFLECT = [
    r"fish$", r"ois$", r"deer$", r"pox$", r"[A-Z].*ese$", r"itis$"]

# table A.1, anglicized
const IRREGULAR_ANG = Dict(
    "beef" => "beefs",
    "brother" => "brothers",
    "child" => "children",
    "cow" => "cows",
    "die" => "dice",
    "ephemeris" => "ephemerides",
    "genie" => "genies",
    "graffito" => "graffiti",
    "money" => "moneys",
    "mongoose" => "mongooses",
    "mythos" => "mythoi",
    "octopus" => "octopuses",
    "ox" => "oxen",
    "soliloquy" => "soliloquies",
    "trilby" => "trilbys")

# table A.1, classical
const IRREGULAR_CLS = merge(IRREGULAR_ANG, Dict(
    "beef" => "beeves",
    "brother" => "brethren",
    "cow" => "kine",
    "die" => "dice",
    "money" => "monies",
    "octopus" => "octopodes"))

# table A.2
const WORD_NOINFLECT = Set([
    "bison", "flounder", "pliers", "bream", "gallows", "proceedings",
    "breeches", "graffiti", "rabies", "britches", "headquarters", "salmon",
    "carp", "herpes", "scissors", "chassis", "high-jinks", "sea-bass",
    "clippers", "homework", "series", "cod", "innings", "shears", "contretemps",
    "jackanapes", "species", "corps", "mackerel", "swine", "debris", "measles",
    "moose", "trout", "diabetes", "mews", "tuna", "djinn", "mumps", "whiting",
    "eland", "news", "wildebeest", "elk", "pincers"])

# table A.3 — currently used only for unpluralization
const SINGLE_S = Set([
    "polis", "asbestos", "glottis", "rhinoceros", "canvas", "lens", "chaos",
    "ibis", "marquis", "caddis", "epidermis", "cannabis", "cosmos", "pathos",
    "dais", "bias", "aegis", "digitalis", "gas", "atlas", "ethos", "pelvis",
    "alias", "trellis", "metropolis", "mantis", "bathos", "acropolis",
    "sassafras"])

# table A.5, anglicized
const PRONOUNS_ANG = Dict(
    "I" => "we",
    "me" => "us",
    "myself" => "ourselves",
    "you" => "you",
    "thou" => "you",
    "thee" => "you",
    "yourself" => "yourself",
    "thyself" => "yourself",
    "she" => "they",
    "he" => "they",
    "it" => "they",
    "they" => "they",
    "her" => "them",
    "him" => "them",
# unfortunately we have no way of telling this from the it above
#    "it" => "them",
    "them" => "them",
    "herself" => "themselves",
    "himself" => "themselves",
    "itself" => "themselves",
    "themself" => "themselves",
    "oneself" => "oneselves")

const PRONOUNS_CLS = let
    temp = copy(PRONOUNS_ANG)
    temp["thou"] = temp["thee"] = "ye"
    temp
end

# table A.10 -a to -ae
const A10 = ("alumna", "alga", "vertebra")

# table A.11 -a to -ae classical
const A11 = Set([
    "abscissa", "formula", "medusa", "amoeba", "hydra", "nebula", "antenna",
    "hyperbola", "nova", "aurora", "lacuna", "parabola"])

# table A.12 -a to -as or -ata
const A12 = Set([
    "charisma", "drama", "carcinoma", "trauma", "schema", "oedema", "enema",
    "dogma", "edema", "sarcoma", "bema", "enigma", "magma", "anathema", "gumma",
    "soma", "miasma", "lymphoma", "melisma", "lemma", "diploma", "stoma",
    "stigma"])

# table A.13 -en to -ens | -ena
const A13 = ("stamen", "foramen", "lumen")

# table A.14 -ex to -ices
const A14 = ("codex", "murex", "silex")

# table A.15 -ex to -exes (ang.) or -ices (cls.)
const A15 = Set([
    "apex", "latex", "vertex", "cortex", "pontifex", "vortex", "index",
    "simplex"])

# table A.16 -is to -ises (ang.) or -ides (cls.)
const A16 = ("iris", "clitoris")

# table A.17 -o to -os
const A17 = Set([
    "albino", "generalissimo", "manifesto", "archipelago", "ghetto", "medico",
    "armadillo", "guano", "octavo", "commando", "inferno", "photo", "ditto",
    "jumbo", "pro", "dynamo", "lingo", "quarto", "embryo", "lumbago", "rhino",
    "fiasco", "magneto", "stylo", "macro"])

# table A.18 -o to -os (ang.) or -i (cls)
const A18 = Set([
    "alto", "contralto", "soprano", "basso", "crescendo", "tempo", "canto",
    "solo"])

# table A.19 -on to -a
const A19 = Set([
    "aphelion", "hyperbaton", "perihelion", "asyndeton", "noumenon",
    "phenomenon", "criterion", "organon", "prolegomenon"])

# table A.20 -um to -a
const A20 = Set([
    "agendum", "datum", "extremum", "bacterium", "desideratum", "stratum",
    "candelabrum", "erratum", "ovum"])

# table A.21 -um to -ums | -a
const A21 = Set([
    "aquarium", "interregnum", "quantum", "compendium", "lustrum", "rostrum",
    "consortium", "maximum", "spectrum", "cranium", "medium", "speculum",
    "curriculum", "memorandum", "stadium", "dictum", "millenium", "trapezium",
    "emporium", "minimum", "ultimatum", "enconium", "momentum", "vacuum",
    "gymnasium", "optimum", "velum", "honorarium", "phylum"])

# table A.22 -us to -uses | -i
const A22 = Set([
    "cactus", "focus", "nimbus", "succubus", "fungus", "nucleolus", "torus",
    "genius", "radius", "umbilicus", "incubus", "stylus", "uterus"])

# table A.23 -us to -uses | -us
const A23 = Set([
    "apparatus", "impetus", "prospectus", "cantus", "nexus", "sinus", "coitus",
    "plexus", "status", "hiatus"])

# in the paper cited, A.24. is used for ifrit to ifriti
# a Google ngrams search reveals no use of ifriti therefore we decide to omit
# A.24 even in classical mode

# table A.25 - to -im
const A25 = ("cherub", "goy", "seraph")

# table A.26 -general to -general
const A26 = ("Adjutant", "Lieutenant", "Quartermaster", "Brigadier", "Major")

const SUFFIX_IRREGULAR = Dict(
    r"man$" => (3, "men"),
    r"[lm]ouse$" => (4, "ice"),
    r"tooth$" => (5, "teeth"),
    r"goose$" => (5, "geese"),
    r"foot$" => (4, "feet"),
    r"zoon$" => (4, "zoa"),
    r"[csx]is$" => (2, "es"))

const SUFFIX_CLASSICAL = Dict(
    r"trix$" => (4, "trices"),
    r"eau$" => (3, "eaux"),
    r"ieu$" => (3, "ieux"),
    r"..[iay]nx$" => (2, "nges"))

const SUFFIX_COMMON = Dict(
    r"[cs]h$" => (1, "hes"),
    r"[xs]$" => (0, "es"),
    r"iz$" => (2, "izzes"),
    r"[^i]z$" => (1, "zes"),
    r"[aeiou]y$" => (1, "ys"),
    r"[A-Z].*y$" => (1, "ys"),
    r"y$" => (1, "ies"),
    r"[aeo]lf$" => (1, "ves"),
    r"[^d]eaf$" => (1, "ves"),
    r"arf$" => (1, "ves"),
    r"[nlw]ife" => (2, "ves"))

function fromend(word, suffixlen)
    offset = length(word) - suffixlen
    offset == 0 ? 0 : nextind(word, 0, offset)
end

stem(word, suffixlen) = word[1:fromend(word, suffixlen)]
stem(word, suffixlen, rest) = stem(word, suffixlen) * rest

function suffix_inflect(table, word)
    for (sfx, (offset, newsfx)) in table
        if occursin(sfx, word)
            return stem(word, offset, newsfx)
        end
    end
    nothing
end

function isnoninflecting(word, classical)
    wordmatches(re) = occursin(re, word)
    any(wordmatches, SUFFIX_NOINFLECT) ||
        word ∈ WORD_NOINFLECT ||
        classical && word ∈ A23
end

"""
    pluralize(word; classical=true)

Pluralize a singular noun `word` (given in canonical capitalization) using
heuristics and lists of exceptions. If `word` is not a singular noun, this
function may give strange results.

If `classical` is set to `true`, then the classical (i.e. inherited from Latin
or Greek) pluralization is chosen instead of the anglicized pluralization. As
an example, the classical plural of `"vertex"` is `"vertices"`, but the
anglicized plural is `"vertexes"`. By default, the classical pluralization is
used.

```jldoctest
julia> using EnglishText

julia> pluralize("fox")
"foxes"

julia> pluralize("radius")
"radii"

julia> pluralize("radius", classical=false)
"radiuses"
```
"""
function pluralize(word::String; classical=true)
    words = split(word)
    if isempty(words)
        throw(ArgumentError("Cannot pluralize the empty string."))
    elseif length(words) == 1
        @assert !isempty(words[1])
        pluralize_single_word(words[1], classical)
    else
        lastword = words[end]
        @assert !isempty(lastword)
        firstwords = join(words[1:end-1], " ")
        if lastword == "General" && !(firstwords ∈ A26)
            string(pluralize(firstwords), " ", lastword)
        else
            string(firstwords, " ", pluralize_single_word(lastword, classical))
        end
    end
end
pluralize(word; classical=true) = pluralize(String(word); classical=classical)

pluralize_single_word(word, classical::Bool) =
    pluralize_single_word(String(word), classical)
function pluralize_single_word(word::String, classical::Bool)::String
    if isnoninflecting(word, classical)
        # words that do not inflect
        return word
    elseif haskey(PRONOUNS_CLS, word)
        # pronouns; note that preposition pronoun not handled currently
        return (classical ? PRONOUNS_CLS : PRONOUNS_ANG)[word]
    elseif haskey(IRREGULAR_CLS, word)
        # irregular plurals
        return (classical ? IRREGULAR_CLS : IRREGULAR_ANG)[word]
    end

    # irregular suffix
    irr = suffix_inflect(SUFFIX_IRREGULAR, word)
    if irr !== nothing
        return coalesce(irr)
    end

    # -in-law
    if endswith(word, "-in-law")
        return string(pluralize_single_word(word[1:end-7], classical),
                      "-in-law")
    end

    # classical inflections, retained
    if word ∈ A14 || classical && word ∈ A15
        return stem(word, 2, "ices")
    elseif word ∈ A20 || classical && word ∈ A21
        return stem(word, 2, "a")
    elseif word ∈ A19
        return stem(word, 2, "a")
    elseif word ∈ A10 || classical && word ∈ A11
        return stem(word, 1, "ae")
    end

    # classical inflections, anglicized
    if classical
        cls = suffix_inflect(SUFFIX_CLASSICAL, word)
        if cls !== nothing
            return coalesce(cls)
        end

        if word ∈ A13
            return stem(word, 2, "ina")
        elseif word ∈ A12
            return stem(word, 1, "ata")
        elseif word ∈ A16
            return stem(word, 2, "ides")
        elseif word ∈ A22
            return stem(word, 2, "i")
        elseif word ∈ A18  # nb: A23 handled as non-inflecting
            return stem(word, 1, "i")
        elseif word ∈ A25
            return word * "im"
        end
    end

    # suffixes ending in -ss, -es, -f, -fe, -y, -ys
    com = suffix_inflect(SUFFIX_COMMON, word)
    if com !== nothing
        return coalesce(com)
    end

    # suffixes ending in -o
    if word ∈ A17 || word ∈ A18
        return word * "s"
    elseif occursin(r"[aeiou]o$", word)
        return word * "s"
    elseif occursin(r"o$", word)
        return word * "es"
    end

    # just add s
    word * "s"
end


text:
using CSV, Distributions, HypothesisTests

data1 = CSV.read("../data/machine1.csv", header=false)[:,1]
data2 = CSV.read("../data/machine2.csv", header=false)[:,1]
xBar1, s1, n1 = mean(data1), std(data1), length(data1)
xBar2, s2, n2 = mean(data2), std(data2), length(data2)
delta0 = 0

sP = sqrt( ( (n1-1)*s1^2 + (n2-1)*s2^2 ) / (n1 + n2 - 2) )
testStatistic = ( xBar1-xBar2 - delta0 ) / ( sP * sqrt( 1/n1 + 1/n2) )
pVal = 2*ccdf(TDist(n1+n2 -2), abs(testStatistic))

println("Manually calculated test statistic: ", testStatistic)
println("Manually calculated p-value: ", pVal, "\n")
println(EqualVarianceTTest(data1, data2, delta0))


text:
#=

  Identify the change point in the coal miners data.
  https://pymc-devs.github.io/pymc/tutorial.html
  """
  Consider the following dataset, which is a time series of recorded
  coal mining disasters in the UK from 1851 to 1962
  [R.G. Jarrett. A note on the intervals between coal mining disasters. Biometrika, 66:191–193, 1979.]
  """

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function change_point(y)
    len = length(y)
    alpha = mean(y) # 1/0.581151832460733  # 1/mean(data)
    lambda_1 ~ Exponential(alpha)
    lambda_2 ~ Exponential(alpha)
    tau ~ DiscreteUniform(1,len)

    # Data is according to Poisson with two different lambdas,
    # one before the change point (tau) and one after the change point
    for i in 1:len
        if tau > i
            y[i] ~ Poisson(lambda_1)
        else
            y[i] ~ Poisson(lambda_2)
        end
    end

end

# Coal miners data from https://pymc-devs.github.io/pymc/tutorial.html
data = [4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6, 3, 3, 5, 4, 5, 3, 1,
        4, 4, 1, 5, 5, 3, 4, 2, 5, 2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3,
        0, 0, 1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1, 0, 1, 0, 1, 0,
        0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2, 3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2,
        0, 0, 1, 4, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1]

model = change_point(data)

#=
println("Prior distribution")
chns_prior = sample(model, Prior(), MCMCThreads(), 1000, num_chns)
display(chns_prior)
display(plot(chns_prior))
=#

println("\nPosterior distribution")
num_chns = 2

# 3.17s tau 40.32+/- 3.77 (rhat ~ 1.1!) 10_000 samples
# chns = sample(model, MH(), MCMCThreads(), 10_000, num_chns)
# 3.17s tau 41.5249+/- 3.03 (rhat ~ 1.04!) 30_000 samples
# chns = sample(model, MH(), MCMCThreads(), 30_000, num_chns)

# Too long...
# chns = sample(model, Gibbs(MH(:tau),PG(20,:lamba_1,:lambda_2)), MCMCThreads(), 10000, num_chns)

# 1.2s tau 57.6+/- 32.1(!) (rhat ok) 1_000 samples
# 3.3s tau 55.98+/- 32.0(!) (rhat ok) 10_000 samples
# chns = sample(model, IS(), MCMCThreads(), 10_000, num_chns)

# 16.2s tau 42.2 +/- 8.5 (rhat 1.05..) 1000 samples
# chns = sample(model, PG(20), MCMCThreads(), 1000, num_chns)

# 2.37s tau 41.1+/- 2.62 (rhat > 2!)  with samples: 1_000
# 17.3s tai 41.28+/- 2.84 (rhat ok)  with samples: 10_000
# chns = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chns)

# 33.1s tau 40.87+/- 2.91 (rhat ~ 1.03)
# chns = sample(model, SMC(), MCMCThreads(), 10_000, num_chns)

# # LoadError: MethodError: Cannot `convert` an object of type Float64 to an object of type Symbol
# chns = sample(model, Gibbs(NUTS(0.65,5,:lambda_1,:lambda_2), PG(20,:tau) ), 1000)

# 19.34s tau 41.4+/-5.18 (rhat ok)
chns = sample(model, Gibbs(HMC(0.1,5,:lambda_1,:lambda_2), PG(20,:tau) ), MCMCThreads(), 1000, num_chns)

# Error
# chns = sample(model, Gibbs(HMCDA(0.15,0.65,0.3,:lambda_1,:lambda_2), PG(20,:tau) ), MCMCThreads(), 1000, num_chns) # Error...

display(chns)
# plot(chns)

show_var_dist_pct(chns, :tau)
show_var_dist_pct(chns, :lambda_1,20)
show_var_dist_pct(chns, :lambda_2,20)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
module ModelSelection

    using DataFrames

    include("structs/modelselection_data.jl")
    include("datatypes/modelselection_result.jl")
    include("const.jl")
    include("strings.jl")
    include("utils.jl")
    include("core.jl")

    include("Preprocessing/Preprocessing.jl")
    include("FeatureExtraction/FeatureExtraction.jl")
    include("PreliminarySelection/PreliminarySelection.jl")
    include("AllSubsetRegression/AllSubsetRegression.jl")
    include("CrossValidation/CrossValidation.jl")
    include("Output/Output.jl")

    using ..Preprocessing
    using ..FeatureExtraction
    using ..PreliminarySelection
    using ..AllSubsetRegression
    using ..CrossValidation
    using ..Output

    export ModelSelectionData, ModelSelectionResult, gsr

    export Preprocessing, FeatureExtraction, PreliminarySelection, Output, CrossValidation

end


text:
@static if isdefined(Turing, :CmdStan)
  function DualAveraging(spl::Sampler{<:AdaptiveHamiltonian}, adapt_conf::CmdStan.Adapt, ϵ::Real)
      # Hyper parameters for dual averaging
      γ = adapt_conf.gamma
      t_0 = adapt_conf.t0
      κ = adapt_conf.kappa
      δ = adapt_conf.delta
      return DualAveraging(γ, t_0, κ, δ, DAState(ϵ))
  end
end

@static if isdefined(Turing, :CmdStan)
    function get_threephase_params(adapt_conf::CmdStan.Adapt)
        init_buffer = adapt_conf.init_buffer
        term_buffer = adapt_conf.term_buffer
        window_size = adapt_conf.window
        next_window = init_buffer + window_size - 1
        return init_buffer, term_buffer, window_size, next_window
    end
end


text:
##############################################################################################################################################
# Implements the Glicko algorithm from Mark Glickman as outlined in the papers: http://www.glicko.net/glicko/glicko.pdf (short version)
# and http://www.glicko.net/research/glicko.pdf (technical version).
# In both cases, the original papers use an artificial point scale for comparability with the Elo system, which complicates the calaculations
# considerably. In the below, we are using natural ratings, which can be related back to the Elo scale by a simple scaling factor.
#
# Natural = Elo_style * ln(10) / 400
#
# In addition, Elo are traditionally started at 1000, where in this case, we set zero as the default rating (so negative scores are possible.)
# This is in keeping with the rest of this package, including the Elo.
##############################################################################################################################################

struct Glicko
    c::Float64
    default_rating::Tuple{Float64, Float64, Int64}
    ratings::Dict
end

function Glicko(;c = 0.1, default_rating = (0, 2.0, -1), ratings = Dict{Int64, Tuple}())
    #Set up a new Elo struct, either using key words or defaults if not provided
    return Glicko(c, default_rating, ratings)
end

#Fit the Glicko rating to some new data
function fit!(glicko::Glicko, original_games::DataFrame)

    games = original_games[!, (1:5)]

    #Duplicate whole rating dataframe with reversed results
    games_duped = dupe_for_rating(games)

    #Split into days and deal with all players in them separately
    for day_games in groupby(games_duped, :Period)
        #Add glicko predictions to each day
        day_games = DataFrame(day_games)
        day = day_games.Period[1]

        #Groupby Player1 including results
        for player_day_games in groupby(day_games, :P1)
            player = player_day_games.P1[1]
            glicko.ratings[player] = update_player_rating(glicko, day, player, player_day_games[!, 2], player_day_games[!, 3], player_day_games[!, 4])
        end
    end
end

function RD_drift(RD, c, defaultRD, day, last_rated)
    #Return the drifted value of RD given the time the rating  was last updated
    if last_rated < 0
        return defaultRD
    else
        return min(sqrt(RD^2 + ((c^2) * abs(day - last_rated))), defaultRD)
    end
end

#g function from Glicko paper
function g(RD::Float64)
    return 1 / sqrt(1 + (0.304 * RD^2))
end

#E function from Glicko paper
function E(r::Float64, rj::Float64, RDj::Float64)
    return 1 / (1 + exp(g(RDj) * (rj - r)))
end

#Reciprocal of d2 function from Glicko paper
function recipd2(r, rj, RDj)
    return sum(g.(RDj).^2 .* E.(r, rj, RDj) .* (1 .- E.(r, rj, RDj)))
end

function update_player_rating(glicko::Glicko, day::Int64, player::Int64, opponents, player_wins, opponent_wins)

    #Collect a tuple of the player being updated...
    (r, RD, last_rated) = get(glicko.ratings, player, glicko.default_rating)
    #Adjusts current RD based on when last rated
    RD = RD_drift(RD, glicko.c, glicko.default_rating[2], day, last_rated)

    #And a vector of the Tuples of the ratings and RDs for their opponents
    rj, RDj, last_rated_j = collect(zip(get.(Ref(glicko.ratings), opponents, Ref(glicko.default_rating))...))
    #Again, giving drift
    RDj = RD_drift.(RDj, Ref(glicko.c), Ref(glicko.default_rating[2]), Ref(day), last_rated_j)

    r_updated = r + ((1 / (1 / RD^2 + recipd2(r, rj, RDj))) * sum(g.(RDj) .* (player_wins .- ((player_wins + opponent_wins) .* E.(r, rj, RDj)))))
    RD_updated = sqrt(1/((1 / RD^2) + recipd2(r, rj, RDj)))

    return (r_updated, RD_updated, day)
end

#Returns the predicted result for any two players in the existing Glicko rating dictionary
function predict(m::Glicko, i::Integer, j::Integer; rating_day::Integer = 0, kwargs...)
    r, RD = rating(m, i, rating_day = rating_day)
    rj, RDj = rating(m, j, rating_day = rating_day)
    return 1 / (1 + exp(g(sqrt(RDj^2 + RD^2)) * (rj - r)))
end

#Returns the rating of a player on a day. If the rating day is before the last day rated, we just take their RD from that day
function rating(m::Glicko, i::Int64; rating_day::Integer = 0)
    r, RD, last_rated = get(m.ratings, i, m.default_rating)

    if rating_day > last_rated
        RD = RD_drift(RD, m.c, m.default_rating[2], rating_day, last_rated)
    end

    return r, RD
end

function rating(m::Glicko, i::Int64, j::Int64; rating_day::Integer = 0)
    return rating(m, i, rating_day = rating_day), rating(m, j, rating_day = rating_day)
end

function display_ratings(m::Glicko; rating_day::Integer = 0)
    a = collect(keys(m.ratings))
    b = map.(x -> glicko.ratings[x][1], a)
    c = map.(x -> glicko.ratings[x][2], a)
    d = map.(x -> glicko.ratings[x][3], a)

    if rating_type == "elo"
        b = convert_natural_to_elo.(b)
        c = convert_natural_to_elo.(c)
    end

    players = DataFrame(player_code = a, player_rating = b, player_sd = c, last_day_played = d)
    return sort!(players, :player_rating, rev = true)
end


text:
using DICOMClient

using FileIO
using ImageMagick
using Images
using JSON3
using ReferenceTests
using Test

@testset "Basic test" begin
    study = "1.2.156.112536.1.2116.222245005117096205.13585270080.1"
    series = "1.2.156.112536.1.2116.222245005117096205.13585270080.2"
    instance = "1.2.156.112536.1.2116.222245005117096205.13585271850.5"
    frames = 1:1
    anonymous_auth = DICOMClient.AnonymousAuth()
    all_auths = [
        anonymous_auth,
    ]
    for auth in all_auths
        base_url = DICOMClient.BaseURL("https://dicomcloud.azurewebsites.net/api")
        client = DICOMClient.Client(base_url, auth)
        @test DICOMClient.get_auth(client) == auth
        @test DICOMClient.get_base_url(client) == base_url
        all_studies_json = DICOMClient.get_all_studies(client)
        @test all_studies_json isa AbstractVector
        all_series_json = DICOMClient.get_all_series(client;
                                                     study = study)
        @test all_series_json isa AbstractVector
        all_instances_json = DICOMClient.get_all_instances(client;
                                                           study = study,
                                                           series = series)
        @test all_instances_json isa AbstractVector
        single_instance_metadata = DICOMClient.get_single_instance_metadata(client;
                                                                            study = study,
                                                                            series = series,
                                                                            instance = instance)
        @test all_instances_json isa AbstractVector
        images = DICOMClient.get_frame_images(Val(:jpeg),
                                              client;
                                              study = study,
                                              series = series,
                                              instance = instance,
                                              frames = frames)
        julia_version_to_expected_hashes = Dict(
            v"1.5" => [0x0f638ae086d9a3c8, 0x0f638ae086d9a3c8],
            v"1.6" => [0x0f638ae086d9a3c8, 0x0f638ae086d9a3c8],
            v"1" =>   [0x0f638ae086d9a3c8, 0x0f638ae086d9a3c8],
        )
        expected_hashes = _get_expected_hashes(julia_version_to_expected_hashes, Base.VERSION)
        @test images isa AbstractVector
        @test length(images) == 2
        for j = 1:2
            @test images[j] isa Matrix{<:RGB}
            @test size(images[j]) == (910, 1260)
            @test hash(images[j]) == expected_hashes[j]
        end
        Base.shred!(auth)
        Base.shred!(client)
    end
    for i in 1:length(all_auths)
        Base.shred!(all_auths[i])
    end
end


text:
module Echidna

using Serialization, Statistics, Random, LinearAlgebra

include("workflow.jl")
include("simulate.jl")
include("compare.jl")
end # module


text:
using Random
using GroupsCore

struct DirectPower{Gr,N,GEl<:GroupsCore.GroupElement} <: GroupsCore.Group
    group::Gr

    function DirectPower{N}(G::GroupsCore.Group) where {N}
        @assert N > 1
        return new{typeof(G),N,eltype(G)}(G)
    end
end

struct DirectPowerElement{GEl,N,Gr<:GroupsCore.Group} <: GroupsCore.GroupElement
    elts::NTuple{N,GEl}
    parent::DirectPower{Gr,N,GEl}
end

DirectPowerElement(
    elts::AbstractVector{<:GroupsCore.GroupElement},
    G::DirectPower,
) = DirectPowerElement(ntuple(i -> elts[i], _nfold(G)), G)

_nfold(::DirectPower{Gr,N}) where {Gr,N} = N

Base.one(G::DirectPower) =
    DirectPowerElement(ntuple(_ -> one(G.group), _nfold(G)), G)

Base.eltype(::Type{<:DirectPower{Gr,N,GEl}}) where {Gr,N,GEl} =
    DirectPowerElement{GEl,N,Gr}

function Base.iterate(G::DirectPower)
    itr = Iterators.ProductIterator(ntuple(i -> G.group, _nfold(G)))
    res = iterate(itr)
    @assert res !== nothing
    elt = DirectPowerElement(first(res), G)
    return elt, (iterator = itr, state = last(res))
end

function Base.iterate(G::DirectPower, state)
    itr, st = state.iterator, state.state
    res = iterate(itr, st)
    res === nothing && return nothing
    elt = DirectPowerElement(first(res), G)
    return elt, (iterator = itr, state = last(res))
end

function Base.IteratorSize(::Type{<:DirectPower{Gr,N}}) where {Gr,N}
    Base.IteratorSize(Gr) isa Base.HasLength && return Base.HasShape{N}()
    Base.IteratorSize(Gr) isa Base.HasShape && return Base.HasShape{N}()
    return Base.IteratorSize(Gr)
end

Base.size(G::DirectPower) = ntuple(_ -> length(G.group), _nfold(G))

GroupsCore.order(::Type{I}, G::DirectPower) where {I<:Integer} =
    convert(I, order(I, G.group)^_nfold(G))

GroupsCore.ngens(G::DirectPower) = _nfold(G)*ngens(G.group)

function GroupsCore.gens(G::DirectPower, i::Integer)
    k = ngens(G.group)
    ci = CartesianIndices((k, _nfold(G)))
    @boundscheck checkbounds(ci, i)
    r, c = Tuple(ci[i])
    tup = ntuple(j -> j == c ? gens(G.group, r) : one(G.group), _nfold(G))
    return DirectPowerElement(tup, G)
end

function GroupsCore.gens(G::DirectPower)
    N = _nfold(G)
    S = gens(G.group)
    tups = [ntuple(j->(i == j ? s : one(s)), N) for i in 1:N for s in S]

    return [DirectPowerElement(elts, G) for elts in tups]
end

Base.isfinite(G::DirectPower) = isfinite(G.group)

function Base.rand(
    rng::Random.AbstractRNG,
    rs::Random.SamplerTrivial{<:DirectPower},
)
    G = rs[]
    return DirectPowerElement(rand(rng, G.group, _nfold(G)), G)
end

GroupsCore.parent(g::DirectPowerElement) = g.parent
GroupsCore.parent_type(::Type{DirectPowerElement{GEl,N,Gr}}) where {GEl,N,Gr} =
    DirectPower{Gr,N,GEl}

Base.:(==)(g::DirectPowerElement, h::DirectPowerElement) =
    (parent(g) === parent(h) && g.elts == h.elts)

Base.hash(g::DirectPowerElement, h::UInt) = hash(g.elts, hash(parent(g), h))

Base.deepcopy_internal(g::DirectPowerElement, stackdict::IdDict) =
    DirectPowerElement(Base.deepcopy_internal(g.elts, stackdict), parent(g))

Base.inv(g::DirectPowerElement) = DirectPowerElement(inv.(g.elts), parent(g))

function Base.:(*)(g::DirectPowerElement, h::DirectPowerElement)
    @assert parent(g) === parent(h)
    return DirectPowerElement(g.elts .* h.elts, parent(g))
end

GroupsCore.order(::Type{I}, g::DirectPowerElement) where {I<:Integer} =
    convert(I, reduce(lcm, (order(I, h) for h in g.elts), init = one(I)))

Base.isone(g::DirectPowerElement) = all(isone, g.elts)

Base.show(io::IO, G::DirectPower) =
    print(io, "Direct $(_nfold(G))-th power of $(G.group)")
Base.show(io::IO, g::DirectPowerElement) =
    print(io, "( ", join(g.elts, ", "), " )")


text:
######################################################################
# General purpose macros
# -----
# Licensed under MIT License

export @await
macro await(exprs...)
    esc(:(fetch(@sync @async begin; $(exprs...); end)))
end

export @sym_str
macro sym_str(str)
    esc(:(Symbol($str)))
end

# Like curry, but applies it to all first-level method calls in a block.
export @curry
"""`@curry exprs...`

 Injects `exprs` as the first arguments in every first-level method calls of the last argument. Like [`curry`](@ref),
 but does not generate an anonymous method instance that can be passed to other methods.

 # Example

 ```julia
 logid = 42
 @curry "Log " logid ": " begin
    println(420)
    println(69)
 end
 # prints "Log 42: 420"
 # prints "Log 42: 69"
 ```
 """
macro curry(curry...)
    if length(curry) < 2
        throw(ArgumentError("no curry-expressions"))
    end

    curry = collect(curry)
    exec = pop!(curry)

    posargs, kwargs = curry_splitargs(curry)

    if exec isa Expr
        if exec.head === :block
            for expr ∈ exec.args
                if expr isa Expr
                    if is_call_expr(expr)
                        curry_call!(expr, posargs, kwargs)
                    elseif is_call_assignment(expr)
                        curry_call!(expr.args[2], posargs, kwargs)
                    end
                end
            end
        elseif is_call_expr(exec)
            curry_call!(exec, posargs, kwargs)
        elseif is_call_assignment(exec)
            curry_call!(exec.args[2], posargs, kwargs)
        end
    end

    esc(exec)
end

function curry_splitargs(args)
    posargs = []
    kwargs  = []
    for arg ∈ args
        if curry_iskwarg(arg)
            push!(kwargs, arg)
        else
            push!(posargs, arg)
        end
    end
    posargs, kwargs
end
curry_iskwarg(_) = false
curry_iskwarg(arg::Expr) = (arg.head === :(=)) || (arg.head === :call && arg.args[1] === :(=>))

function curry_call!(call::Expr, pos, kw)
    @assert call.head === :call
    splice!(call.args, 2:1, pos)
    append!(call.args, curry_makekwarg.(kw))
    call
end
function curry_makekwarg(arg::Expr)
    if arg.head === :(=)
        Expr(:kw, arg.args[1], arg.args[2])
    elseif arg.head === :call && arg.args[1] === :(=>)
        Expr(:kw, arg.args[2], arg.args[3])
    else
        throw(ArgumentError("Unknown expression $arg"))
    end
end

export @with
"""`@with resources... block`
 Produces code like:
 ```
 let resources...
    try
        block
    finally
        close.(resources)
    end
 end
 ```

 This is useful for resources implementing `Base.close` but no callbacked resource opener.
 """
macro with(exprs::Union{Symbol, Expr}...)
    block = last(exprs)
    assignments = with_assignments(__source__, exprs[1:lastindex(exprs)-1])
    vars = extract_with_varname.(assignments)

    closers = Expr(:., :close, Expr(:tuple, Expr(:tuple, vars...)))
    wrap    = Expr(:try, block, false, false, closers)
    esc(Expr(:let, Expr(:block, assignments...), wrap))
end

with_assignments(__source__, exprs) = (with_assignment(__source__, expr, Ref{Int}(0)) for expr in exprs)
function with_assignment(__source__, expr::Expr, counter::Ref{<:Integer})
    if expr.head === :(=)
        expr
    else
        counter[] += 1
        Expr(:(=), Symbol("$(__source__.file):$(__source__.line):$(counter[])"), expr)
    end
end
with_assignment(_, var::Symbol, _) = Expr(:(=), var, var)

function extract_with_varname(expr::Expr)
    @assert expr.head === :(=) "Not an assignment"
    expr.args[1]
end

export @once
macro once(cond::Union{Bool, Expr}, expr)
    label = QuoteNode(Symbol("$(__source__.file):$(__source__.line)"))
    esc(quote
        if $label ∉ KirUtil.ONCE_CONDITIONS && $cond
            push!(KirUtil.ONCE_CONDITIONS, $label)
            $expr
        end
    end)
end

const ONCE_CONDITIONS = Set{Symbol}()


text:
@test mdinfo(:symmetric | :posdef) != nothing

@test mdinfo(isloaded) != nothing


text:
struct FakeNaturalDiffWithIsApprox  # For testing overloading isapprox(::Tangent) works:
    x
end
function Base.isapprox(c::Tangent, d::FakeNaturalDiffWithIsApprox; kwargs...)
    return isapprox(c.x, d.x, kwargs...)
end
function Base.isapprox(d::FakeNaturalDiffWithIsApprox, c::Tangent; kwargs...)
    return isapprox(c.x, d.x, kwargs...)
end

@testset "check_result.jl" begin
    @testset "_test_add!!_behaviour" begin
        check = ChainRulesTestUtils._test_add!!_behaviour

        check(10.0, 2.0)
        check(11.0, ZeroTangent())
        check([10.0, 20.0], @thunk([2.0, 0.0]))

        check(12.0, InplaceableThunk(X̄ -> error("Should not have in-placed"), @thunk(2.0)))

        check([10.0, 20.0], InplaceableThunk(X̄ -> (X̄[1] += 2.0; X̄), @thunk([2.0, 0.0])))
        @test fails() do
            check([10.0, 20.0], InplaceableThunk(X̄ -> (X̄[1] += 3.0; X̄), @thunk([2.0, 0.0])))
        end
    end

    @testset "test_approx" begin
        @testset "positive cases" begin
            test_approx(1.0, 1.0)
            test_approx(1.0 + im, 1.0 + im)
            test_approx(1.0, 1.0 + 1e-10)  # isapprox _behaviour
            test_approx((1.5, 2.5, 3.5), (1.5, 2.5, 3.5 + 1e-10))

            test_approx(ZeroTangent(), 0.0)

            test_approx([1.0, 2.0], [1.0, 2.0])
            test_approx([[1.0], [2.0]], [[1.0], [2.0]])

            test_approx(@thunk(10 * 0.1 * [[1.0], [2.0]]), [[1.0], [2.0]])

            test_approx(@not_implemented(""), rand(3))
            test_approx(rand(3), @not_implemented(""))
            test_approx(@not_implemented("a"), @not_implemented("a"))

            test_approx(
                Tangent{Tuple{Float64,Float64}}(1.0, 2.0),
                Tangent{Tuple{Float64,Float64}}(1.0, 2.0),
            )

            diag_eg = Diagonal(randn(5))
            test_approx( # Structual == Structural
                Tangent{typeof(diag_eg)}(; diag=diag_eg.diag),
                Tangent{typeof(diag_eg)}(; diag=diag_eg.diag),
            )
            test_approx( # Structural == Natural
                Tangent{typeof(diag_eg)}(; diag=diag_eg.diag),
                diag_eg,
            )

            T = (a=1.0, b=2.0)
            test_approx(
                Tangent{typeof(T)}(; a=1.0), Tangent{typeof(T)}(; a=1.0, b=ZeroTangent())
            )
            test_approx(
                Tangent{typeof(T)}(; a=1.0),
                Tangent{typeof(T)}(; a=1.0 + 1e-10, b=ZeroTangent()),
            )

            test_approx(
                Tangent{FakeNaturalDiffWithIsApprox}(; x=1.4),
                FakeNaturalDiffWithIsApprox(1.4),
            )
            test_approx(
                FakeNaturalDiffWithIsApprox(1.4),
                Tangent{FakeNaturalDiffWithIsApprox}(; x=1.4),
            )
        end
        @testset "negative case" begin
            @test fails(() -> test_approx(1.0, 2.0))
            @test fails(() -> test_approx(1.0 + im, 1.0 - im))
            @test fails(() -> test_approx((1.5, 2.5, 3.5), (1.5, 2.5, 4.5)))

            @test fails(() -> test_approx(ZeroTangent(), 20.0))
            @test fails(() -> test_approx(10.0, ZeroTangent()))

            @test fails(() -> test_approx([1.0, 2.0], [1.0, 3.9]))
            @test fails(() -> test_approx([[1.0], [2.0]], [[1.1], [2.0]]))

            @test fails(() -> test_approx(@thunk(10 * [[1.0], [2.0]]), [[1.0], [2.0]]))

            @test fails(() -> test_approx(@not_implemented("a"), @not_implemented("b")))
        end
        @testset "type negative" begin
            @test fails() do  # these have different primals so should not be equal
                test_approx(
                    Tangent{Tuple{Float32,Float32}}(1.0f0, 2.0f0),
                    Tangent{Tuple{Float64,Float64}}(1.0, 2.0),
                )
            end
            @test fails() do
                test_approx((1.0, 2.0), Tangent{Tuple{Float64,Float64}}(1.0, 2.0))
            end
        end

        @testset "TestIterator" begin
            data = randn(3)
            iter1 = TestIterator(data, Base.HasLength(), Base.HasEltype())
            iter2 = TestIterator(data, Base.HasLength(), Base.EltypeUnknown())
            test_approx(iter2, iter1)

            iter3 = TestIterator(data .+ 1e-10, Base.HasLength(), Base.HasEltype())
            test_approx(iter3, iter1)

            iter_bad = TestIterator(data .+ 010, Base.HasLength(), Base.HasEltype())
            @test fails(() -> test_approx(iter_bad, iter1))
        end
    end
end


text:
# Copyright (c) 2020 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

function get_loggrad!(grada::AbstractVector{T},
                      gradA::AbstractArray{T},
                      grady::AbstractArray{T},
                      betat0::Vector{T}, betat1::Vector{T},
                      Nt::Integer,
                      a::AbstractVector{T},
                      A::AbstractMatrix{T},
                      y::AbstractArray{T},
                      alpha::Matrix{T},
                      logML::T
                     ) where {T <: AbstractFloat}
  Ns = size(A,1)
  C = logML*Nt + log(T(Nt))
  for t = Nt-1:-1:1
    for j = 1:Ns
      for k = 1:Ns
        ajk = A[j,k]
        if ajk != -T(Inf)
          betaA = ajk + betat1[k] + y[t+1,k]
          betat0[j] = logadd(betat0[j], betaA)
          betaA += alpha[t,j]
          grady[t+1,k] = logadd(grady[t+1,k], betaA)
          gradA[j,k]   = logadd(gradA[j,k]  , betaA)
        end
      end
    end
    betat0, betat1 = betat1, betat0
    fill!(betat0,-T(Inf))
  end
  @views grady[1,:] .=  betat1 .+ alpha[1,:]
  grady .= -exp.(grady.-C)
  gradA .= -exp.(gradA.-C)
  @views grada .= grady[1,:]
  return grada, gradA, grady
end

# with constrained path
function get_loggrad!(grada,
                      gradA,
                      grady::AbstractArray{T},
                      betat0::Vector{T}, betat1::Vector{T},
                      Nt::Integer,
                      t2IJ::Vector{Pair{Vector{D},Vector{D}}},
                      A::AbstractMatrix{T},
                      y::AbstractArray{T},
                      alpha::Matrix{T},
                      logML::T) where {D<:Integer,
                                       T <: AbstractFloat}
  C = logML*Nt+log(T(Nt))
  Ns = size(A,1)
  for t = Nt-1:-1:1
    I,J = t2IJ[t]
    for k in eachindex(I)
      i,j = I[k],J[k]
      betaA = A[i,j] + betat1[j] + y[t+1,j]
      betat0[i] = logadd(betat0[i], betaA)
      betaA += alpha[t,i]
      grady[t+1,j] = logadd(grady[t+1,j], betaA)
      #gradA[i,j]   = logadd(gradA[i,j]  , betaA)
    end
    betat0, betat1 = betat1, betat0
    fill!(betat0,-T(Inf))
  end
  idx0 = unique(t2IJ[1][1])
  for i in idx0
    grady[1,i] = betat1[i] + alpha[1,i]
  end
  grady .= -exp.(grady.-C)
  #gradA .= -exp.(gradA.-C)
  return grada, gradA, grady
end

function get_loggrad(Nt::Integer,
                     a,
                     A::AbstractMatrix{T},
                     y::AbstractArray{T},
                     alpha::Matrix{T},
                     logML::T
                    ) where {T <: AbstractFloat}
  Ns = size(A,1)
  grada = typeof(a) <: Array{T} ? -T(Inf) .* ones(T,Ns) : NoTangent()
  gradA = typeof(a) <: Array{T} ? -T(Inf) .* ones(T,Ns,Ns) : NoTangent()
  grady = -T(Inf) .* ones(T,size(y,1),Ns)
  betat1 = zeros(T,Ns)
  betat0 = -T(Inf).*ones(T,Ns)
  get_loggrad!(grada,gradA,grady,betat0,betat1,Nt,a,A,y,alpha,logML)
end

# y is a tensor Nt2 × Ns × Nb
function get_loggrad(
                     Nts::Vector{D},
                     as::AbstractMatrix{T},
                     As::AbstractArray{T,3},
                     ys::AbstractArray{T,3},
                     alphas::Vector{Matrix{T}},
                     logMLs::Vector{T}
                    ) where {T <: AbstractFloat, D <: Integer}
  Ns = size(As,1)
  Nb = size(ys,3)
  grada = -T(Inf) .* ones(T,Ns,Nb)
  gradA = -T(Inf) .* ones(T,Ns,Ns,Nb)
  grady = -T(Inf) .* ones(T,size(ys,1),Ns,Nb)
  betat1 = zeros(T,Ns)
  betat0 = -T(Inf).*ones(T,Ns)
  for i = 1:Nb
    get_loggrad!(
                 view(grada,:,i),
                 view(gradA,:,:,i),
                 view(grady,:,:,i),
                 betat0,betat1,
                 Nts[i],
                 view(as,:,i),
                 view(As,:,:,i),
                 view(ys,:,:,i),
                 alphas[i],logMLs[i])
    fill!(betat1, zero(T))
    fill!(betat0, -T(Inf))
  end
  return grada,gradA,grady
end

# y is a tensor Nt2 × Ns × Nb, constrained path
function get_loggrad(
                     Nts::Vector{D},
                     as::Vector{K},
                     A::AbstractMatrix{T},
                     y::AbstractArray{T,3},
                     alphas::Vector{Matrix{T}},
                     logMLs::Vector{T}
                    ) where {T <: AbstractFloat, Z, K<:AbstractVector{Z}, D <: Integer}
  Ns = size(A,1)
  Nb = size(y,3)
  grada = NoTangent()
  gradA = NoTangent()
  grady = -T(Inf) .* ones(T,size(y,1),Ns,Nb)
  betat1 = zeros(T,Ns)
  betat0 = -T(Inf).*ones(T,Ns)
  for i = 1:Nb
    get_loggrad!(
                 grada,
                 gradA,
                 view(grady,:,:,i),
                 betat0,betat1,
                 Nts[i],as[i],
                 A,
                 view(y,:,:,i),
                 alphas[i],logMLs[i])
    fill!(betat1, zero(T))
    fill!(betat0, -T(Inf))
  end
  return grada,gradA,grady
end


text:
module StaleA

stale(x) = rand(1:8)
stale(x::Int) = length(digits(x))

not_stale(x::String) = first(x)

use_stale(c) = stale(c[1]) + not_stale("hello")
build_stale(x) = use_stale(Any[x])

# force precompilation
build_stale(37)
stale('c')

end


text:
using CoinbasePro
using Test
using DataFrames

@testset "CoinbasePro" begin

for file in filter(x->occursin("test_", x), readdir("."))
  include(file)
end

end


text:
# test basic forward convolve, see IIF issue #477

# using Revise

using Test
using IncrementalInference
using Statistics

@testset "Test basic convolution result..." begin

function forwardConvolve(X0::Array{Float64,2}, model)

  fg = initfg()

  addVariable!(fg, :x0, ContinuousScalar)
  initManual!(fg, :x0, X0)
  addVariable!(fg, :x1, ContinuousScalar)

  addFactor!(fg, [:x0;:x1], model)

  ## TODO -- dont use name here, add API to just use z2 here
  return approxConv(fg, :x0x1f1, :x1)
end

## Start

# first numerical values -- samples from the marginal of X0
z1 = Normal(0,0.1)
X0 = rand(z1, 1,100)

## predict -- project / conv
# 0 -> 1 seconds
# make approx function
z2 = Normal(11,1.0) # odo
statemodel = LinearRelative( z2 )
X1_ = forwardConvolve(X0, statemodel)

## measure -- product of beliefs, using `ApproxManifoldProducts.jl`

predX1 = manikde!(X1_, ContinuousScalar)
z3 = Normal(9.5,0.75)
measX1 = manikde!(reshape(rand(z3,100),1,:), ContinuousScalar)

# do actual product
posterioriX1 = predX1 * measX1
X1 = getPoints(posterioriX1)

## predict, 1->2 seconds
z4 = Normal(8,2.0) # odo
statemodel = LinearRelative( z4 )
X2_ = forwardConvolve(X1, statemodel)

@test size(X2_) == (1,100)
@test 15 < Statistics.mean(X2_) < 25

end


text:
@discord_object struct EmbedImage
    url::String
    proxy_url::String
    height::Int
    width::Int
end


text:
# NOTATER: Det er øpv å skrive følgande (gitt at TemporalType er "using"-introdusert i moder-module.
#   NyStruct <: TemporalType
# TODO: Du kan kjøre T<:TemporalType !

mutable struct Tile_map{T} <: Representation
    _range
    _number_of_intervals::Unsigned
    _all_SAT #::Vector{T}
    function Tile_map{T}(interval::Union{Nothing, Tuple{<:Real, <:Real}};
                         N::Union{Int,Missing}=missing,    the_sat=missing) where {T}

        if !ismissing(the_sat) && !isnothing(interval)
            throw(ArgumentError("Cannot have a single SAT for a Euclidean interval"))
        end
        !isnothing(interval) && ismissing(N)   &&   throw(ArgumentError("Interv is supplied but not N"))
        # Trur denne er unødvendig, sjå proxy-ctor under.
        #
        the_sat_vector = []
        if isnothing(interval)
            if ismissing(the_sat)
                the_sat_vector = [T()]
            else
                the_sat_vector = [the_sat]
            end
            N = 1
        else
            # Dersom N is nothing, kast en error.
            for i in first(interval):1/N:last(interval)
                push!(the_sat_vector, T(i))
            end
            # Remove last?
            pop!(the_sat_vector)
            # the_map = Dict(i => T(10*i) for i=range(first(interval)), last(interval), length=1/N)
        end

        new{T}(interval, N, the_sat_vector)
    end
end

function Tile_map{T}(;the_sat=missing) where {T}
    if ismissing(the_sat)
        the_sat = T()
    end
    Tile_map{T}(nothing, the_sat=the_sat)
end

function dimention_of_NRES(arg::Tile_map)
    if isnothing(arg._range)
        return 0
    else
        return 1
    end
    # (plan: kan den være rekursiv for composite NRES? )
end

"""
    function map_to_SAT(nres::Tile_map, coordinate::Real)

Map scalar coordinate to NRES-til set.
"""
function map_to_SAT(nres::Tile_map; coordinate::Union{Real, Missing}=missing)
    # TODO TODO  Kkjlør TDD
    #
    # kjør full TDD her, dene må testas!
    #
    #
    #
    if isnothing(nres._range)
        return nres._all_SAT[1]
    end
    ismissing(coordinate) && throw(ArgumentError("map_to_SAT(1D_nres) with no coord? Error!"))

    # TODO TEST_coordinate is missing when range is not missing: Should raise error!

    # Assuming 1D coordinate (c ∈ ℜ¹) within range of NRES:
    if coordinate < first(nres._range) || coordinate > last(nres._range)
        throw(ArgumentError(string("Coordinate ", coordinate, " is not in range ", nres._range)))
    end

    # 0: create a linspace with the same steps
    value_list = range(first(nres._range), last(nres._range), length=nres._number_of_intervals)
    # 1: find index of first entry above coordinate value
    indeks = findfirst(value_list.>=coordinate)
    # 2: Use this index when conllecting corresponding SAT from NRES map:
    nres._all_SAT[indeks]
end

export dimention_of_NRES

max_stars_repo_name: leikanger/NRES

max_issues_repo_name: leikanger/NRES

max_forks_repo_name: leikanger/NRES


text:
using DiffEqBase: bisection
using Test

# basic tests
# result of bisection should never be zero,
# should have sign same as left limit
# and nexfloat in tdir should be zero or same sign of right limit
r = bisection(sin, (3.0, 4.0), 1.0)
@test !iszero(sin(r))
@test sin(r) > 0
@test sin(nextfloat(r)) <= 0

r = bisection(sin, (4.0, 3.0), -1.0)
@test !iszero(sin(r))
@test sin(r) < 0
@test sin(prevfloat(r)) >= 0

# this example has a big zero interval, makes it go inside the inner loop
for i in 1:10
	f = (x) -> (10.0^(-i))*(x^3)
	r = bisection(f, (-1.0, 1.0), 1.0)
	@test !iszero(f(r))
	@test f(r) < 0
	@test f(nextfloat(r)) >= 0

	r = bisection(f, (1.0, -1.0), -1.0)
	@test !iszero(f(r))
	@test f(r) > 0
	@test f(prevfloat(r)) <= 0
end


text:
rgbs = IsoLumc
luvs = Luv.(IsoLumc)
lchs = LCHuv.(IsoLumc)

using Test
ls = comp1.(luvs)
@test all(comp1.(luvs) .== comp1.(lchs))

scatter(comp2.(luvs),comp3.(luvs),color=rgbs,markersize=9)
scatter(hue.(lchs),comp2.(lchs),color=rgbs,markersize=9,proj=:polar)

max_stars_repo_name: VLAB-ExperimentSystem/VLab

max_issues_repo_name: VLAB-ExperimentSystem/VLab

max_forks_repo_name: VLAB-ExperimentSystem/VLab


text:
function vtk_xml_write_header(vtk::DatasetFile)
    xroot = create_root(vtk.xdoc, "VTKFile")
    set_attribute(xroot, "type", vtk.grid_type)
    set_attribute(xroot, "version", "1.0")
    if IS_LITTLE_ENDIAN
        set_attribute(xroot, "byte_order", "LittleEndian")
    else
        set_attribute(xroot, "byte_order", "BigEndian")
    end
    set_attribute(xroot, "header_type", string(HeaderType))
    if vtk.compression_level > 0
        set_attribute(xroot, "compressor", "vtkZLibDataCompressor")
    end
    xroot::XMLElement
end

# Returns the "extent" attribute required for structured (including rectilinear)
# grids.
extent_attribute(Ni, Nj, Nk, ::Nothing=nothing) =
    "0 $(Ni - 1) 0 $(Nj - 1) 0 $(Nk - 1)"

function extent_attribute(Ni, Nj, Nk, extent::Array{T}) where T <: Integer
    length(extent) == 6 || throw(ArgumentError("extent must have length 6."))
    (extent[2] - extent[1] + 1 == Ni) &&
    (extent[4] - extent[3] + 1 == Nj) &&
    (extent[6] - extent[5] + 1 == Nk) ||
    throw(ArgumentError("extent is not consistent with dataset dimensions."))
    join(extent, " ")
end

# Number of cells in structured grids.
# In 3D, all cells are hexahedrons (i.e. VTK_HEXAHEDRON), and the number of
# cells is (Ni-1)*(Nj-1)*(Nk-1). In 2D, they are quadrilaterals (VTK_QUAD), and in
# 1D they are line segments (VTK_LINE).
function num_cells_structured(Ni, Nj, Nk)
    Ncls = one(Ni)
    for N in (Ni, Nj, Nk)
        Ncls *= max(1, N - 1)
    end
    Ncls
end

# Accept passing coordinates as a tuple.
vtk_grid(filename, coords::NTuple; kwargs...) =
    vtk_grid(filename, coords...; kwargs...)