text:
using MomentClosure
using MomentClosure: define_M, define_μ
using ModelingToolkit: value
using SymbolicUtils: expand
using Test
using Catalyst

@parameters t, c₁, c₂, c₃, c₄, Ω
@variables X(t), Y(t)

S_mat = [ 1 -1  1 -1;
         -1  1  0  0]
a = [c₁*X*Y*(X-1)/Ω^2, c₂*X, c₃*Ω, c₄*X]
rn = ReactionSystemMod(t, [X, Y], [c₁, c₂, c₃, c₄, Ω], a, S_mat)

μ = define_μ(2,4)
M = define_M(2,4)
sys = generate_central_moment_eqs(rn, 2, 4)
expr1 = sys.odes.eqs[1].rhs

closed_eqs = moment_closure(sys, "zero")
@test closed_eqs.closure[M[2,2]] == 0
expr1 = closed_eqs.odes.eqs[1].rhs
expr2 = c₃*Ω + M[1,1]*c₁*μ[1,0]*(Ω^-2) + M[1,1]*c₁*(Ω^-2)*(μ[1,0]- 1) + c₁*M[2,0]*μ[0,1]*(Ω^-2) +
        c₁*μ[0,1]*μ[1,0]*(Ω^-2)*(μ[1,0] - 1) - c₂*μ[1,0] - c₄*μ[1,0]
expr2 = simplify(value.(expr2))
@test isequal(expand(expr1), expand(expr2))

# check that deterministic_IC is working with central moments
ic_values = Dict(deterministic_IC([2, 5], closed_eqs))
@test length(ic_values) == 5
@test ic_values[μ[1,0]] == 2 && ic_values[μ[0,1]] == 5 && ic_values[M[1,1]] == 0

closed_eqs = moment_closure(sys, "normal")
@test isequal(closed_eqs.closure[M[0,4]], 3*M[0,2]^2)

closed_eqs= moment_closure(sys, "log-normal")
expr1 = closed_eqs.closure[M[1,2]]
expr2 = μ[1,0]*μ[0,1]^2*(1.0+M[0,2]*μ[0,1]^-2)*(1.0 + M[1,1]*(μ[0,1]^-1)*(μ[1,0]^-1))^2 -
        M[0,2]*μ[1,0] - μ[1,0]*μ[0,1]^2 - 2*M[1,1]*μ[0,1]
@test isequal(expr1, expr2)

closed_eqs = moment_closure(sys, "poisson")
@test isequal(closed_eqs.closure[M[3,0]], μ[1,0])

closed_eqs = moment_closure(sys, "gamma")
expr1 = closed_eqs.closure[M[2,1]]
expr2 = M[2,0]*μ[0,1] + μ[0,1]*μ[1,0]^2 + 2*M[1,1]*μ[1,0] + 2*M[1,1]*M[2,0]*μ[1,0]^-1 -
    M[2,0]*μ[0,1] - μ[0,1]*μ[1,0]^2 - 2*M[1,1]*μ[1,0]
@test isequal(expr1, expr2)
expr1

closed_eqs = moment_closure(sys, "derivative matching")
expr1 = closed_eqs.closure[sys.M[0,4]]
expr2 = μ[0,1]^4*(M[0,2]+μ[0,1]^2)^-6*(M[0,3]+μ[0,1]^3+3*M[0,2]*μ[0,1])^4 - μ[0,1]^4 -
    6*M[0,2]*μ[0,1]^2 - 4*M[0,3]*μ[0,1]
@test isequal(expr1, expr2)

sys = generate_raw_moment_eqs(rn, 2)
closed_eqs = moment_closure(sys, "zero")
expr1 = closed_eqs.closure[μ[3,0]]
expr2 = -2*μ[1,0]^3 + 3*μ[1,0]*μ[2,0]
@test isequal(expr1, expr2)

closed_eqs = moment_closure(sys, "normal")
expr1 = closed_eqs.closure[sys.μ[1,2]]
expr2 = 2*μ[0,1]*μ[1,1] + μ[0,2]*μ[1,0] - 2*μ[1,0]*μ[0,1]^2
@test isequal(expr1, expr2)
expr1 = closed_eqs.odes.eqs[5].rhs
expr2 = c₂*μ[1,0] + 2*c₂*μ[1,1] + c₁*μ[0,1]*μ[2,0]/Ω^2 - c₁*μ[1,1]/Ω^2 - 4*c₁*μ[1,1]^2/Ω^2 + 4*c₁*μ[0,1]*μ[1,1]/Ω^2 -
    2*c₁*μ[0,1]*μ[1,0]^2/Ω^2 + 2*c₁*μ[0,2]*μ[1,0]/Ω^2 -2*c₁*μ[0,2]*μ[2,0]/Ω^2 + 2*c₁*μ[1,0]*μ[1,1]/Ω^2 -
    4*c₁*μ[1,0]*μ[0,1]^2/Ω^2 + 4*c₁*μ[0,1]^2*μ[1,0]^2/Ω^2
expr2 = simplify(value.(expr2))
@test isequal(expand(expr1), expr2)

# check that deterministic_IC is working with raw moments
ic_values = Dict(deterministic_IC([2, 5], closed_eqs))
@test length(ic_values) == 5
@test ic_values[μ[1,0]] == 2 && ic_values[μ[0,1]] == 5 && ic_values[μ[1,1]] == 10

closed_eqs = moment_closure(sys, "log-normal")
expr1 = closed_eqs.closure[μ[1,3]]
expr2 = μ[0,1]^-6 * μ[0,2]^3 * μ[1,0]^-2 * μ[1,1]^3
@test isequal(expr1, expr2)

closed_eqs = moment_closure(sys, "poisson")
expr1 = closed_eqs.closure[μ[4,0]]
expr2 = μ[1,0] + 6 * μ[1,0]^4 + 3*μ[2,0]^2 + 4*μ[1,0]*μ[3,0] -12*μ[2,0]*μ[1,0]^2
@test isequal(expr1, expr2)

closed_eqs = moment_closure(sys, "gamma")
expr1 = closed_eqs.closure[μ[0,3]]
expr2 = μ[0,1]^3 + 3*μ[0,1]*(μ[0,2]-μ[0,1]^2) + 2*μ[0,1]^-1*(μ[0,2]-μ[0,1]^2)^2
@test isequal(expr1, expr2)

closed_eqs = moment_closure(sys, "derivative matching")
expr1 = closed_eqs.closure[sys.μ[0,3]]
expr2 = μ[0,1]^-3*μ[0,2]^3
@test isequal(expr1, expr2)


text:
#using SegwayRobot

function dynmodel!(x::AbstractVector, u::AbstractVector, dt::Number, nn_params)
    # Segway definition
    Mb = 1.76 # kg # Masse des Körpers
    Mw = 0.147 # kg # Masse der Räder
    R = 0.07 # m #Radius der Räder
    cz = 0.077 # m #Schwerpunkt des Körpers Höhe über Radachse
    b = 0.1985/2 # m #Half the distance between the wheels (OO_w in paper)
    Ixx = (0.166^2+0.21^2)*Mb/12 + Mb*cz^2 # kg m^2 #Inertia of the body
    Iyy = (0.072^2+0.21^2)*Mb/12 + Mb*cz^2 # kg m^2 #Inertia of the body
    Izz = (0.072^2+0.166^2)*Mb/12 # kg m^2 #Inertia of the body
    Iwa = Mw*R^2/2 # kg m^2 #Inertia of the wheels about their axis
    Iwd = Mw*b^2 # kg m^2 #Inertia of the wheels about their axis
    seg = Segway_unmutable(Mb, Mw, R, cz, b, Ixx, Iyy, Izz, Iwa, Iwd)

    # Evaluate Neural Net for the torque
    tau = eval_nn_segway(nn_params, vcat(x[4:7], u))

    # Evaluate segway dynamics
    segway_rk4!(x, tau, dt, seg)

    return x
end

@inline function eval_nn_segway(w, x)
    for i=1:2:length(w)-2
        x = tanh.(w[i] * x .+ w[i+1])
    end
    return w[end-1] * x .+ w[end]
end

"""
State vector is [x0, y0, phi, alpha, dalpha, v, dphi, ir, il]
"""
function segway_cost(xt::AbstractVector, u::AbstractVector)
    cost = (xt[4]^2 + 0.1*xt[5]^2
            + 10*xt[2]^2 + 10*xt[1]^2 # or squared here...
            + 0.1* xt[7]^2
            + 0.01*(u'u)) # punish for falling over

    reward = 0.0 #reward for moving along the x-axis

    if xt[4]^2 > 1 #check if lager than pi/2
        terminal = 1
        cost *= 1.2
    else
        terminal = 0
    end

    return cost-reward, terminal
end

ext: jl

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
using Gridap
using GridapEmbedded
using Test

# Define BC functions
println("Defining Boundary conditions")
H = 0.41
Um= 1.0
u_in(x, t) = VectorValue(1.5 * Um * x[2] * (H - x[2]) / ((H / 2)^2), 0.0)
u_noSlip(x, t) = VectorValue(0.0, 0.0)
u_in(t::Real) = x -> u_in(x, t)
u_noSlip(t::Real) = x -> u_noSlip(x, t)
u_noSlip(x::VectorValue) = u_noSlip(x,0.0)

f(x) = VectorValue(0, 0)
g(x) = 0.0

# R = 0.3
# pmin = Point(0,0)
# pmax = Point(1,1)
# n = 20
# partition = (n,n)

# geo1 = disk(R,x0=Point(0.5,0.5))
# geo2 = ! geo1

# bgmodel = CartesianDiscreteModel(pmin,pmax,partition)
# dp = pmax - pmin

## Build geometry
# Cilinder
R = 0.05
C = Point(0.2,0.2)
cylinder = disk(R,x0=C)
# Flag
C1 = Point(0.2,0.19)
d1 = VectorValue(0.4,0.0)
d2 = VectorValue(0.0,0.02)
C2 = Point(0.6,0.21)
n1 = VectorValue(1.0,0.0)
n2 = VectorValue(0.0,1.0)
p1 = plane(x0=C1,v=-n1)
p2 = plane(x0=C1,v=-n2)
p3 = plane(x0=C2,v=n1)
p4 = plane(x0=C2,v=n2)
flag = intersect(intersect(p1,p2),intersect(p3,p4))
not_fluid = union(cylinder,flag)
fluid_geo = !(not_fluid,name="fluid")

# Background model
pmin = Point(0.0,0.0)
pmax = Point(2.2,0.41)
partition = (120,60)
function map(x)
  γ=3.0
  if x[2]<=0.2
    f0 = exp.(γ*0.2/0.41)
    y = (f0 - exp.(-γ*(x[2]-0.2)/0.41)) / (f0-1.0) * 0.2
  else
    f0 = exp.(γ*0.21/0.41)
    y = 0.2 + (exp.(γ*(x[2]-0.2)/0.41)-1.0) / (f0-1.0) * 0.21
  end
  return VectorValue(x[1],y)
end
bgmodel = CartesianDiscreteModel(pmin,pmax,partition,map=map)

# Boundary labels
labels = get_face_labeling(bgmodel)
add_tag_from_tags!(labels,"wall",[2,4,5,6])
add_tag_from_tags!(labels,"inlet",[1,3,7])
add_tag_from_tags!(labels,"outlet",[8])

dp = pmax - pmin
const h = dp[1]/80

#cutgeo = cut(bgmodel,geo2)
cutgeo = cut(bgmodel,fluid_geo)
model = DiscreteModel(cutgeo)

strategy = AggregateAllCutCells()
aggregates = aggregate(strategy,cutgeo)

trian = Triangulation(bgmodel)
trian_Ω = Triangulation(cutgeo)
trian_Γd = EmbeddedBoundary(cutgeo)
trian_Γn = BoundaryTriangulation(model,"outlet")

order = 2
quad_Ω = CellQuadrature(trian_Ω,2*order)
quad_Γd = CellQuadrature(trian_Γd,2*order)
quad_Γn = CellQuadrature(trian_Γn,2*order)

n_Γd = get_normal_vector(trian_Γd)
n_Γn = get_normal_vector(trian_Γn)

Vstd = TestFESpace(
  reffe=:QLagrangian,
  conformity=:H1,
  valuetype=VectorValue{2,Float64},
  model=model,
  order=order,
  dof_space=:physical,
  dirichlet_tags=["inlet", "wall"])

Vser = TestFESpace(
  reffe=:SLagrangian,
  conformity=:L2, # we don't neet to impose continuity since we only use the cell dof basis / shapefuns
  valuetype=VectorValue{2,Float64},
  model=model,
  order=order,
  dof_space=:physical)

Qstd = TestFESpace(
  reffe=:PLagrangian,
  conformity=:L2,
  valuetype=Float64,
  model=model,
  order=order-1,
  dof_space=:physical)

V = AgFEMSpace(Vstd,aggregates,Vser)
Q = AgFEMSpace(Qstd,aggregates)

U = TrialFESpace(V,[u_in(0.0),u_noSlip(0.0)])
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V,Q])
X = MultiFieldFESpace([U,P])

function A_Ω(x,y)
  u, p = x
  v, q = y
  ∇(v)⊙∇(u) - q*(∇⋅u) - (∇⋅v)*p
end

function B_Ω(y)
  v, q = y
  v⋅f - q*g
end

const γ = order*(order+1)

function A_Γd(x,y)
  u, p = x
  v, q = y
  (γ/h)*v⋅u - v⋅(n_Γd⋅∇(u)) - (n_Γd⋅∇(v))⋅u + (p*n_Γd)⋅v + (q*n_Γd)⋅u
end

function B_Γd(y)
  v, q = y
  (γ/h)*v⋅u_noSlip - (n_Γd⋅∇(v))⋅u_noSlip + (q*n_Γd)⋅u_noSlip
end

function B_Γn(y)
  v, q = y
  v⋅(n_Γn⋅∇(u)) - (n_Γn⋅v)*p
end

t_Ω = AffineFETerm(A_Ω,B_Ω,trian_Ω,quad_Ω)
t_Γd = AffineFETerm(A_Γd,B_Γd,trian_Γd,quad_Γd)
t_Γn = FESource(B_Γn,trian_Γn,quad_Γn)
op = AffineFEOperator(X,Y,t_Ω,t_Γd)#,t_Γn)
uh, ph = solve(op)

uh_Ω = restrict(uh,trian_Ω)
ph_Ω = restrict(ph,trian_Ω)

writevtk(trian_Ω,"stokes.vtu",cellfields=["vh"=>uh_Ω,"ph"=>ph_Ω])

#colors = color_aggregates(aggregates,bgmodel)
#writevtk(trian,"trian",celldata=["aggregate"=>aggregates,"color"=>colors],cellfields=["uh"=>uh,"ph"=>ph])
#writevtk(model,"model")
#writevtk(trian_Ω,"trian_O",cellfields=["uh"=>uh_Ω,"ph"=>ph_Ω])
#writevtk(trian_Γd,"trian_Gd",cellfields=["normal"=>n_Γd])
#writevtk(trian_Γn,"trian_Gn",cellfields=["normal"=>n_Γn])

#eu = u - uh_Ω
#ep = p - ph_Ω

#l2(u) = u⊙u
#h1(u) = ∇(u)⊙∇(u) + l2(u)

# eul2 = sqrt(sum( integrate(l2(eu),trian_Ω,quad_Ω) ))
# euh1 = sqrt(sum( integrate(h1(eu),trian_Ω,quad_Ω) ))
# epl2 = sqrt(sum( integrate(l2(ep),trian_Ω,quad_Ω) ))

# @test eul2 < 1.e-6
# @test euh1 < 1.e-6
# @test epl2 < 1.e-6


text:
# Autogenerated wrapper script for Libgpg_error_jll for armv7l-linux-musleabihf
export libgpg_error

JLLWrappers.@generate_wrapper_header("Libgpg_error")
JLLWrappers.@declare_library_product(libgpg_error, "libgpg-error.so.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgpg_error,
        "lib/libgpg-error.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Code used to generate robust control problems.
# Save in Testproblems/RobustControl with an accurante
# solution from MOSEK

# using JuMP, Mosek, Base.Test, HelperFunctions, JLD
# Create number of Lyapunov stable controller and a number of Hinf performance controller problems
nnLyapunov = 0
nnHinf = 25

function LyapunovStabilityLMI(A,Bu)
  n = size(A,1)
  nu = size(Bu,2)
  AP  = kron(eye(n),A)
  PAt = kron(A,eye(n))
  BuY = kron(eye(n),Bu)
  YtBut = kron(Bu,eye(n))
  # turn YtBut*vec(Y') into YtBut* vec(Y) by permuting cols of YtBut
  YtBut = HelperFunctions.transposeVectorized(YtBut,n,nu) #provide dimension of Y'

  # to get strict inequality and prevent trivial solution P=0
  ba = vec(-1e-8*speye(n))

  return [AP+PAt BuY+YtBut], ba
end

function HinfPerformanceLMI(A,Bu,Bw,Cz,Dzw,Dzu)
  n,nw = size(Bw)
  nz,nu = size(Dzu)
  # -----------------------------------
  # LMI): Hinf performance and stability
  # -----------------------------------

  # at first create the subblocks of the LMI left hand side
  # [A1     0       A2
  #  0      A3       0
  # A4      0       A5]
  # ---- A1 -----
  AP  = kron(eye(n),A)
  PAt = kron(A,eye(n))
  BuY = kron(eye(n),Bu)
  YtBut = kron(Bu,eye(n))
  YtBut = transposeVectorized(YtBut,n,nu) #provide dimension of Y'
  A1 = [AP+PAt BuY+YtBut zeros(n^2,1)] #will be multiplied by [vec(P),vec(Y),γ]^T
  # ---- A2 -----
  PCzt = kron(Cz,eye(n))
  YtDzut = kron(Dzu,eye(n))
  YtDzut = transposeVectorized(YtDzut,n,nu) #provide dimension of Y'
  A2 = [PCzt YtDzut zeros(n*nz,1)]
  # ---- A3 -----
  A3 = [zeros(nw^2,n^2+nu*n) -vec(eye(nw))]
  # ---- A4 -----
  CzP = kron(eye(n),Cz)
  DzuY = kron(eye(n),Dzu)
  A4 = [CzP DzuY zeros(nz*n,1)]
  # ---- A5 -----
  A5 = [zeros(nz^2,n^2+nu*n) -vec(eye(nz))]

  # Assemble big LMI block, s.t. Aa * [vec(P),vec(Y),γ]^T = vec(LMI)
  Aa = [A1[1:n,:];zeros(nw,size(A1,2));A4[1:nz,:]]
  for i = 1:n-1
    Aa = [Aa;A1[i*n+1:(i+1)*n,:];zeros(nw,size(A1,2));A4[i*nz+1:(i+1)*nz,:]]
  end
  for i = 0:nw-1
    Aa = [Aa;zeros(n,size(A1,2));A3[i*nw+1:(i+1)*nw,:];zeros(nz,size(A1,2))]
  end
  for i = 0:nz-1
    Aa = [Aa;A2[i*n+1:(i+1)*n,:];zeros(nw,size(A1,2));A5[i*nz+1:(i+1)*nz,:]]
  end

  Ba = [zeros(n,n) -Bw zeros(n,nz);
      -Bw' zeros(nw,nw) -Dzw';
      zeros(nz,n) -Dzw zeros(nz,nz)]
  # add matrix with small eigenvalues to right hand side to achieve a strict inequality
  Ba = Ba - (1e-8)*eye(size(Ba,1))
  ba = vec(Ba)
  # create right hand side matrix B, all constant entries of LMI
  return Aa,ba
end

rng = MersenneTwister(125)
dirPath = "../DataFiles/Julia/RobustControl/"
!ispath(dirPath) && mkdir(dirPath)

iii = 1
 while iii <=(nnLyapunov+nnHinf)

  # create random LTI system
  # choose number of states x(t)
  n = rand(rng,3:6)
  # choose number of outputs for z(t)
  nz = rand(rng,3:6)
  #choose number of inputs for u(t) and w(t)
  nu = rand(rng,3:6)
  nw = rand(rng,3:6)
  # create system matrices (create system in such a way that a stabilizing controller K exists)
  A = randn(rng,n,n)

  # choose Bu in such a way that (A,Bu) is controllable
  isControllable = false
  Bu = 0
  while !isControllable
    Bu = sprandn(rng,n,nu,0.8)
    # check controllability
    Con = Bu
    for kkk=1:n-1
      res = A^kkk*Bu
      Con = [Con res]
    end
    isControllable = (rank(full(Con)) == n)
  end
  Bw = eye(n,nw)
  Cz = randn(rng,nz,n)
  Dzw = randn(rng,nz,nw)
  Dzu = eye(nz,nu)

  # LMI 1): Lyapunov stability: PA'+AP+Y'Bu'+BuY < 0 (NOTE STRICT INEQUALITY)
  if iii <=nnLyapunov
    Aa, ba = LyapunovStabilityLMI(A,Bu)
    Pa = spzeros(n^2+nu*n,n^2+nu*n)
    qa = spzeros(n^2+nu*n,1)
    nLMI = n
    # add P > 0 constraint
    Aa = [Aa; -speye(n^2) spzeros(n^2,nu*n)]
    ba = [ba;zeros(n^2)]
  else
    # LMI 2): Hinf performance:
    Aa,ba = HinfPerformanceLMI(A,Bu,Bw,Cz,Dzw,Dzu)
    nLMI = n+nw+nz
    Pa = spzeros(n^2+nu*n+1,n^2+nu*n+1)
    qa = [spzeros(n^2+nu*n,1);1]
    # add P > 0 constraint
    Aa = [Aa; -speye(n^2) spzeros(n^2,nu*n+1)]
    ba = [ba;vec(-1e-8*eye(n))]
  end
  r = 0.
    # specify cone
  Kf = 0
  Kl = 0
  Kq = []
  Ks = [nLMI^2 n^2]

  #solve accurately once with mosek
  model = Model(solver=MosekSolver())
  @variable(model, P[1:n,1:n], SDP)
  @variable(model, Y[1:nu,1:n])
  @variable(model, g)
  @objective(model, Min , g)
  LMI = [A*P+P*A'+Bu*Y+Y'*Bu' Bw P*Cz'+Y'*Dzu';
          Bw' -g*eye(nw) Dzw';
          Cz*P+Dzu*Y Dzw -g*eye(nz)]
  @SDconstraint(model, LMI <= -1e-8*eye(nLMI))
  status = JuMP.solve(model)
  objTrue = getobjectivevalue(model)
  P = getvalue(P)
  Y = getvalue(Y)
  g = getvalue(g)
  K = Y/P
  Ac = A+Bu*K
  eig(Ac)

  # compute once gain in other format
  model2 = Model(solver=MosekSolver())
  @variable(model2, S1[1:n,1:n], SDP)
  @variable(model2, S2[1:nLMI,1:nLMI], SDP)
  @variable(model2, P2[1:n,1:n], SDP)
  @variable(model2, Y2[1:nu,1:n])
  if iii <= nnLyapunov
    @objective(model2, Min , 0)
    x2 = [vec(P2);vec(Y2)]
  else
    @variable(model2, g2)
    @objective(model2, Min , g2)
    x2 = [vec(P2);vec(Y2);g2]
  end
  s2 = [vec(S2);vec(S1)]
  @constraint(model2,Aa*x2+s2.==ba)
  status = JuMP.solve(model2)
  objTrue = getobjectivevalue(model2)
  P2 = getvalue(P2)
  Y2 = getvalue(Y2)
  K2 = Y2/P2
  g2 = getvalue(g2)
  Ac2 = A+Bu*K2

  # check eigenvalues of c.l. system A+Bu*F
    unstableK = maximum(map(x->real(x),eig(Ac2)[1])) >= 0.

  # discard all the weird cases with stalling solver or unstable eigenvalues
  if status != :Optimal || unstableK
    warn("Bad problem: Try again!")
    continue
  end

  # Some further tests for development
  # @test norm(Ac-Ac2,Inf) < 1e-5
  # @test abs(g-g2) < 1e-5
  # @test norm(vec(getvalue(LMI))-(Aa*[vec(P);vec(Y);g]-ba),Inf) < 1e-5
  # @test maximum(map(x->real(x),eig(Ac)[1])) <= 0.
  # @test maximum(map(x->real(x),eig(Ac2)[1])) <= 0.
  nr = "$(iii)"
  if iii < 10
    nr = "0$(iii)"
  end
  iii <= nnLyapunov ? problemName = "Lyap"*nr : problemName = "Hinf"*nr
  fn = "RobustControl"*nr*".jld"
  problemType = "LMI Robust Control Problem"
  JLD.save(dirPath*fn,"n",size(Aa,2),"m",size(Aa,1),"A",Aa,"b",ba,"P",Pa,"q",qa,"r",r,"objTrue",objTrue,"K",K2,"problemType",problemType,"problemName",problemName,"Kf",Kf,"Kl",Kl,"Kq",Kq,"Ks",Ks,"nu",nu,"nz",nz,"nw",nw,"nx",n)
  println("$(iii)/$(nnLyapunov+nnHinf) completed!")
  iii+=1
end


text:
using ReducedBasisMethods
using Test

@testset "ReducedBasisMethods.jl" begin
    include("parameter_tests.jl")
    include("parameterspace_tests.jl")
    include("trainingset_tests.jl")
end


text:
@fenicsclass Geometry

#functions necessary for creating meshes from geometrical objects.
#2D objects below
Circle(centre,radius) = Geometry(mshr.Circle(centre,radius))
Rectangle(corner1,corner2)=Geometry(mshr.Rectangle(corner1,corner2))
Ellipse(centre,horizontal_semi_axis,vertical_semi_axis,fragments)=Geometry(mshr.Ellipse(centre,horizontal_semi_axis,vertical_semi_axis,fragments))

#3d objects below
Box(corner1,corner2) = Geometry(mshr.Box(corner1,corner2))
Cone(top,bottom,bottom_radius,slices::Int)=Geometry(mshr.Cone(top,bottom,bottom_radius,slices))
Sphere(centre,radius) = Geometry(mshr.Sphere(centre,radius))

#we generate the mesh of some geometrical objects
generate_mesh(geom_object::Geometry,size::Int)=Mesh(mshr.generate_mesh(geom_object.pyobject,size))
export Circle,Rectangle,Ellipse,Box,Cone,Sphere,generate_mesh

#operator overloading for Geometry types so we can create composite shapes
+(geom_object1::Geometry, geom_object2::Geometry) = Geometry(geom_object1.pyobject.__add__(geom_object2.pyobject))
-(geom_object1::Geometry, geom_object2::Geometry) = Geometry(geom_object1.pyobject.__sub__(geom_object2.pyobject))
*(geom_object1::Geometry, geom_object2::Geometry) = Geometry(geom_object1.pyobject.__mul__(geom_object2.pyobject))

set_subdomain(object::Geometry,degree,domain2) = fenicspycall(object, :set_subdomain,degree,domain2.pyobject)
"""
Extrude2D

Extrudes a 2D geometry to 3D
"""
Extrude2D(object::Geometry,thickness) = Geometry(mshr.Extrude2D(object.pyobject,thickness))
export  set_subdomain, Extrude2D

CSGUnion(geom_object1::Geometry, geom_object2::Geometry) = Geometry(mshr.CSGUnion(geom_object1.pyobject,geom_object2.pyobject))
CSGIntersection(geom_object1::Geometry, geom_object2::Geometry) = Geometry(mshr.CSGIntersection(geom_object1.pyobject,geom_object2.pyobject))
CSGDifference(geom_object1::Geometry, geom_object2::Geometry) = Geometry(mshr.CSGDifference(geom_object1.pyobject,geom_object2.pyobject))

CSGScaling(geom_object1::Geometry, scale) = Geometry(mshr.CSGScaling(geom_object1.pyobject,scale))

CSGRotation(geom_object1::Geometry, angle) = Geometry(mshr.CSGRotation(geom_object1.pyobject,angle))
CSGRotation(geom_object1::Geometry, point, angle) = Geometry(mshr.CSGRotation(geom_object1.pyobject,point, angle))
CSGRotation(geom_object1::Geometry, point_axis, point_center, angle) = Geometry(mshr.CSGRotation(geom_object1.pyobject, point_axis, point_center, angle))

CSGTranslation(geom_object1::Geometry, translation_point) = Geometry(mshr.CSGTranslation(geom_object1.pyobject, translation_point))

export CSGUnion, CSGIntersection, CSGDifference, CSGScaling, CSGRotation, CSGTranslation


text:
"""
    collatz_sequence(n)

Collatz conjecture: start with any positive integer n. The next term is
obtained as follows:
	If n term is even, the next term is: n / 2 .
	If n is odd, the next term is: 3 * n + 1.
The conjecture states the sequence will always reach 1 for any starting value n.
"""

function collatz_sequence(n)
	sequence = [n]
	while n != 1
		n = n % 2 == 1 ? 3 * n + 1 : div(n, 2)
		append!(sequence, n)
	end
	return sequence
end

max_stars_repo_name: ashwani-rathee/Julia

max_issues_repo_name: ashwani-rathee/Julia

max_forks_repo_name: ashwani-rathee/Julia


text:
@testset "Cellsets" begin

  # needs expanding
  cellset_result = Cellset("Cellset 1")
  @test name(cellset_result) == "Cellset 1"

  cellset_2 = JSON.parse("""
  {
    "ID": "Cellset 2"
  }
  """)

  @test name(Cellset(cellset_2)) == "Cellset 2"
  @test Dimension(cellset_2) != cellset_result

end


text:
"""
    F = JopPow(spc, c, a)

where `F` is the power operator (x/c)^a with domain and range given by `spc::JetSpace`, and scalar values c, a.
"""
JopPow(spc::JetSpace{T}, c = 1, a = 2) where {T} = JopNl(dom = spc, rng = spc, f! = JopPow_f!, df! = JopPow_df!, df′! = JopPow_df′!, s = (c=T(c), a=T(a)))
export JopPow

JopPow_f!(d::AbstractArray, m::AbstractArray; c, a) = d .= (m ./ c).^a
JopPow_df!(δd::AbstractArray, δm::AbstractArray; mₒ, c, a) = δd .= ((a/c) .* (mₒ ./ c).^(a-1)) .* δm
JopPow_df′!(δm::AbstractArray, δd::AbstractArray; mₒ, c, a)= δm .= conj((a/c) .* (mₒ ./ c).^(a-1)) .* δd


text:
using Logging
using LinearAlgebra
using SparseArrays

using ArgParse
using DelimitedFiles
using TimerOutputs

using JuMP

# Solvers
import CPLEX
import Gurobi
import Mosek
import MosekTools

# CPLEX cut parameters
# We set all these to 0 to de-activate all cuts
const CPLEX_CUT_PARAMS = [
    "CPXPARAM_MIP_Cuts_BQP",
    "CPXPARAM_MIP_Cuts_Cliques",
    "CPXPARAM_MIP_Cuts_Covers",
    "CPXPARAM_MIP_Cuts_Disjunctive",
    "CPXPARAM_MIP_Cuts_FlowCovers",
    "CPXPARAM_MIP_Cuts_PathCut",
    "CPXPARAM_MIP_Cuts_Gomory",
    "CPXPARAM_MIP_Cuts_GUBCovers",
    "CPXPARAM_MIP_Cuts_Implied",
    "CPXPARAM_MIP_Cuts_LocalImplied",
    "CPXPARAM_MIP_Cuts_LiftProj",
    "CPXPARAM_MIP_Cuts_MIRCut",
    # "CPXPARAM_MIP_Cut_MCFCut",  # For some reason this one raises an error
    "CPXPARAM_MIP_Cuts_RLT",
    "CPXPARAM_MIP_Cuts_ZeroHalfCut"
]

using CLaP

import Base.RefValue

function compute_analytic_center(sf, optimizer)
    model = Model(optimizer)
    set_optimizer_attribute(model, MOI.Silent(), true)
    m, n = size(sf.A)

    # Build continuous relaxation
    x = Vector{VariableRef}(undef, n)

    for (kidx, k) in sf.cones
        d = MOI.dimension(k)
        x[kidx] .= @variable(model, [1:d] in k)
    end

    @constraint(model, sf.A*x .== sf.b)

    @objective(model, Min, 0)

    optimize!(model)

    xref = value.(x)

    return xref
end

function landp_closure(
    finst::String,
    micp_optimizer::MOI.OptimizerWithAttributes,
    cgcp_optimizer::MOI.OptimizerWithAttributes,
    time_limit::Float64, nrm::Symbol, max_rounds::Int;
    verbose::Bool=true,
    max_refine_rounds::Int=50, conic_feas_tol::Float64=1e-1,
    timer=TimerOutput()
)
    # Read model from file
    cbf = MOI.FileFormats.CBF.Model()
    MOI.read_from_file(cbf, finst)

    # Convert to standard form
    @timeit timer "Build SF" sf = build_standard_form(micp_optimizer, cbf, bridge_type=Float64)

    # Compute analytic center
    xref = compute_analytic_center(sf, micp_optimizer)

    # TODO: add pre-processing option
    @timeit timer "Prepross" extract_implied_integer(sf)
    @info "Pre-processing: $(sum(sf.vartypes)) integer variables ($(sum(sf.vartypes_implied)) implied integers)"

    # TODO: add cut validity checker option
    # x_micp = readdlm(joinpath(@__DIR__, "../cplex_micp/res/$(basename(finst)[1:end-4]).sol"))
    # @info "MICP objective value: $(dot(sf.c, x_micp))"

    # Set some parameters
    m, n = size(sf.A)
    micp = sf.model
    MOI.set(micp, MOI.TimeLimitSec(), time_limit)
    MOI.set(micp, MOI.Silent(), !verbose)

    # Book-keeping
    ncalls = 0
    nrefine = 0
    nscuts_tot = 0
    nkcuts_tot = 0
    nnz_kcut_tot = 0
    nnz_scut_tot =
    nrounds = 0

    # Set callback
    x = sf.var_indices
    tstart = 0.0

    function cblandp(cbdata)
        ncalls += 1

        nrounds > max_rounds && return nothing

        # Get fractional point
        x_ = MOI.get(micp, MOI.CallbackVariablePrimal(cbdata), x)
        z_ = dot(sf.c, x_)  # Current lower bound

        # Refine outer approximation before separating cuts
        @timeit timer "Refine OA" if nrefine < max_refine_rounds
            H = CLaP.extract_conic_infeasibilities(x_, sf.cones) ./ 2
            kflag = false
            for ((kidx, k), η) in zip(sf.cones, H)
                isa(k, CLaP.POLYHEDRAL_CONE) && continue
                @assert isa(k, MOI.SecondOrderCone) "Only (R)SOC are supported (is $k)"

                if η <= -conic_feas_tol / 2

                    # Add K* cut
                    λ = copy(x_[kidx])
                    λ[1] = norm(λ[2:end])
                    λ[2:end] .*= -1
                    λ ./= λ[1]

                    MOI.submit(
                        micp,
                        MOI.UserCut(cbdata),
                        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(λ, x[kidx]), 0.0),
                        MOI.GreaterThan(0.0)
                    )
                    kflag = true
                end
            end

            nrefine += kflag
            kflag && return nothing
        end

        nrounds += 1
        nrefine = 0  # reset refine counter

        # Compute cuts
        Kcuts, Scuts = CLaP.lift_and_project(
            x_, sf, cgcp_optimizer,
            normalization=nrm,
            kcut_pre_check=true,
            kcut_post_check=true,
            strengthen_flag=true,
            timer=timer,
            xref=xref
        )

        # Submit cuts
        nnz_kcut = 0
        for (js, α, β) in Kcuts
            nnz_kcut += length(α)
            MOI.submit(
                micp,
                MOI.UserCut(cbdata),
                MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(α, x[js]), 0.0),
                MOI.GreaterThan(β)
            )
        end

        nnz_scut = 0
        for (js, α, β) in Scuts
            nnz_scut += length(α)
            MOI.submit(
                micp,
                MOI.UserCut(cbdata),
                MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(α, x[js]), 0.0),
                MOI.GreaterThan(β)
            )
        end

        # Book-keeping
        nkcuts = length(Kcuts)
        nscuts = length(Scuts)
        t_ = time() - tstart
        @info "Stats for round $nrounds" nkcuts nscuts z_ t_ nnz_kcut nnz_scut

        nkcuts_tot += nkcuts
        nscuts_tot += nscuts
        nnz_scut_tot += nnz_scut
        nnz_kcut_tot += nnz_kcut

        return nothing
    end

    MOI.set(micp, MOI.UserCutCallback(), cblandp)

    # Solve
    tstart = time()
    @timeit timer "MICP" MOI.optimize!(micp)

    # Result log
    @info "User callback was called $(ncalls[]) times" nscuts_tot nnz_scut_tot nkcuts_tot nnz_kcut_tot

    @info MOI.get(micp, MOI.TerminationStatus())
    @info "Final bound: $(MOI.get(micp, MOI.ObjectiveBound()))"

    return micp
end

function parse_commandline(cl_args)
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--MICPSolver"
            help = "MICP solver"
            arg_type = Symbol
            default = :CPLEX
        "--CGCPSolver"
            help = "CGCP solver"
            arg_type = Symbol
            default = :Gurobi
        "--Normalization"
            help = "Normalization condition"
            arg_type = Symbol
            default = :Standard
        "--Rounds"
            help = "Maximum number of cutting plane rounds"
            arg_type = Int
            default = typemax(Int)
        "--TimeLimit"
            help = "Time limit (in seconds)"
            arg_type = Float64
            default = 1e30
        "finst"
            help = "Instance file (in CBF format)"
            required = true
    end

    return parse_args(cl_args, s)
end

function main(CLARGS)
    cl_args = parse_commandline(CLARGS)
    @info("User options", cl_args)

    BLAS.set_num_threads(1)

    # Set solvers
    if cl_args["MICPSolver"] == :CPLEX
        micp_optimizer = MOI.OptimizerWithAttributes(
            CPLEX.Optimizer,
            "CPX_PARAM_PREIND" => 0,
            "CPX_PARAM_THREADS" => 1,
            # disable all cuts
            [cutparam => -1 for cutparam in CPLEX_CUT_PARAMS]...,
            "CPXPARAM_MIP_Limits_CutsFactor" => 1e30,
            "CPXPARAM_MIP_Strategy_HeuristicFreq" => -1,
            "CPXPARAM_MIP_Limits_Nodes" => 0,
            "CPXPARAM_MIP_Strategy_MIQCPStrat" => 2  # OA algorithm
        )
    elseif cl_args["MICPSolver"] == :Gurobi
        GRB_ENV = Gurobi.Env()  # To avoid checking out multiple Gurobi licenses
        micp_optimizer = MOI.OptimizerWithAttributes(
            () -> Gurobi.Optimizer(GRB_ENV),
            "Threads" => 1,
            "Presolve" => 0,
            "Heuristics" => 0,
            "Cuts" => 0,
            "MIQCPMethod" => 1,
            "NodeLimit" => 1,
            "QCPDual" => 0
        )
    else
        error("MICP solver $(cl_args["MICPSolver"]) is not supported.\n Possible values are `CPLEX` and `Gurobi`.")
    end

    if cl_args["CGCPSolver"] == :CPLEX
        cgcp_optimizer = MOI.OptimizerWithAttributes(
            CPLEX.Optimizer,
            "CPX_PARAM_THREADS" => 1,
            "CPX_PARAM_PREIND" => 0,
            "CPX_PARAM_SCRIND" => 0
        )

    elseif cl_args["CGCPSolver"] == :Gurobi
        GRB_ENV = Gurobi.Env()  # To avoid checking out multiple Gurobi licenses
        cgcp_optimizer = MOI.OptimizerWithAttributes(
            () -> Gurobi.Optimizer(GRB_ENV),
            "Threads" => 1,
            "OutputFlag" => 0,
            "Presolve" => 0,
            "QCPDual" => 0
        )

    elseif cl_args["CGCPSolver"] == :Mosek
        cgcp_optimizer = MOI.OptimizerWithAttributes(
            Mosek.Optimizer,
            "NUM_THREADS" => 1,
            "LOG" => 0
        )
    else
        error("CGCP solver $(cl_args["CGCPSolver"]) is not supported.\n Possible values are `CPLEX`, `Gurobi` and `Mosek`.")
    end

    # Warm-up run
    to = TimerOutput()
    @info "Warming up..."
    with_logger(Logging.NullLogger()) do  # this will de-activate logs
        landp_closure(
            cl_args["finst"],
            micp_optimizer, cgcp_optimizer,
            10.0, cl_args["Normalization"], 2,
            verbose=false, timer = to
        )
    end

    # Real run
    to = TimerOutput()
    landp_closure(
        cl_args["finst"],
        micp_optimizer, cgcp_optimizer,
        cl_args["TimeLimit"], cl_args["Normalization"], cl_args["Rounds"],
        timer=to
    )

    print_timer(to)
    println()

    return nothing
end

# Run the main method if script is called directly
if abspath(PROGRAM_FILE) == @__FILE__
    main(ARGS)
end

max_stars_repo_name: mtanneau/CLaP

max_issues_repo_name: mtanneau/CLaP

max_forks_repo_name: mtanneau/CLaP


text:
using PyPlot, Bridge
λ(t) = 6/(t + 2)
T = 10.0
λmax = 3.0

P = InhomogPoisson(λ)
tt = sample(ThinningAlg(λmax), T, P)

plot(Bridge.piecewise(SamplePath(tt, collect(1:length(tt))))..., label="N(t)")
legend()


text:
using Term
import Term.progress: track, ProgressBar, update, start, stop
import Term.consoles: clear
import Term.color: RGBColor

tprint(hLine("progress bars"; style="blue"))

"""
The easiest way to use Term's progress bars is through the
`track` function. This just takes an iterable object as input
when you're doing a for loop over the object. It will print
out a visualization of progress along the for loop:
"""

myvec = collect(1:100)
for i in track(myvec)
    # do stuff
    sleep(.025)
end

"""
`track` creates and updates a `ProgressBar` object for you.
But you can do it manually if you want more control over it
"""
# create progress bar
pbar = ProgressBar(;N=250, description="Manual pbar", width=120)
# loop da loop
for i in 1:250
    # do stuff
    sleep(.007)

    # manually update the progress bar
    update(pbar)
end
stop(pbar)

"""
As you can see, progress bars have parameters to set their appearance.
The same parameters can be used for both ProgressBar and track.
Here's more:
"""

print("\n")
trk(x) = track(x;
    description="[red bold]Choose your colors![/red bold]",
    expand=true,  #  fill the screen
    update_every=5,  # don't update at every iteration
    columns=:detailed,  # print more info
    colors = [
        RGBColor("(.3, .3, 1)"),
        RGBColor("(1, 1, 1)"),
        RGBColor("(.9, .3, .3)"),
    ]
)

print("\n")
for i in trk(1:500)
    sleep(0.005)
end

"""
Each bit of information (the description, the progress bar, the number of
completed iterations...) is a `column`. A progress bar can be made up
of different sets of columns, and you can create your own.

Term provides three different presets to display different ammounts of information:
"""

print("\n")
for (level, color) in zip((:minimal, :default, :detailed), ("yellow1", "green","blue"))
    for i in track(1:100; description="[bold $color italic]$level[/bold $color italic]", columns=level, width=150)
        sleep(.005)
    end
end

"""
As you can see the display is adjusted so that it always has the right
width regardless of how many columns there are. The only thing
left to see is how to create your own column.
"""

import Term.progress: AbstractColumn, DescriptionColumn, BarColumn
import Term.style: apply_style  # to style strings
import Term.measure: Measure  # to define the column's size

# define a new column type
struct MyCol <: AbstractColumn
    segments::Vector
    measure::Measure
    style::String
end

# constructor
MyCol(style::String) = MyCol([], Measure(10, 1), style)

# define a function to update the column when the progress bar is updated
# it must return a string
Term.progress.update(col::MyCol, i::Int, args...)::String = return apply_style("[$(col.style)]$i[/$(col.style)]")

# okay let's create a progress bar with out columns
cols = [
    DescriptionColumn("MY columns!!!"),
    BarColumn(),  # this will show the actual bar
    MyCol("bold red"),
]

for i in track(1:100; columns=cols)
    sleep(.01)
end

"""
Our simple column simply keeps track of the number of iterations in the loop,
but you can go crazy and make all sorts of columns! Check the docs for more info.

Sometimes you want to create a lot of progress bars, but don't want to
clutter the console space too much, make them transient!
"""

print("\n")
for i in 1:3
    for j in track(1:250; description="[yellow bold]Transient pbars![/yellow bold]", transient=true, width=250)
        sleep(0.001)

    end
end
tprintln("[bright_blue bold]poof![/bright_blue bold] [underline bright_blue]They disappeared[/]")

tprint(hLine("Done"; style="green"))


text:
# """
#     strain2x2tto3v!(v::FVec{T}, t::FMat{T}) where {T}

# Convert a matrix of 2x2 strain components  into a 3-component vector.
# """
# function strain2x2tto3v!(v::FVec{T}, t::FMat{T}) where {T}
#     v[1] = t[1,1];
#     v[2] = t[2,2];
#     v[3] = t[1,2] + t[2,1];
#     return v
# end

# """
#     strain3vto2x2t!(t::FMat{T}, v::FVec{T}) where {T}

# Convert a strain 3-vector to a  matrix of 2x2 strain components (symmetric tensor).
# """
# function strain3vto2x2t!(t::FMat{T}, v::FVec{T}) where {T}
#     t[1,1] = v[1];
#     t[2,2] = v[2];
#     t[1,2] = v[3]/2.;
#     t[2,1] = v[3]/2.;
#     return t
# end

# """
#     strain3x3tto6v!(v::FVec{T}, t::FMat{T}) where {T}

# Convert a matrix of 3x3 strain components to a 6-component vector.
# """
# function strain3x3tto6v!(v::FVec{T}, t::FMat{T}) where {T}
#     v[1] = t[1,1];
#     v[2] = t[2,2];
#     v[3] = t[3,3];
#     v[4] = t[1,2] + t[2,1];
#     v[5] = t[1,3] + t[3,1];
#     v[6] = t[3,2] + t[2,3];
#     return v
# end

# """
#     strain6vto3x3t!(t::FMat{T}, v::FVec{T}) where {T}

# Convert a strain 6-vector to a  matrix of 3x3 strain components (symmetric tensor)..
# """
# function strain6vto3x3t!(t::FMat{T}, v::FVec{T}) where {T}
#     t[1,1] = v[1];
#     t[2,2] = v[2];
#     t[3,3] = v[3];
#     t[1,2] = v[4]/2.;
#     t[2,1] = v[4]/2.;
#     t[1,3] = v[5]/2.;
#     t[3,1] = v[5]/2.;
#     t[3,2] = v[6]/2.;
#     t[2,3] = v[6]/2.;
#     return t
# end

# """
#     strain3x3tto9v!(v::FVec{T}, t::FMat{T}) where {T}

# Convert a matrix of 3x3 strain components to a 9-component vector.

# The  strain components are in the order
#      ex, ey, ez, gxy, gyx, gyz, gzy, gxz, gzx
# """
# function strain3x3tto9v!(v::FVec{T}, t::FMat{T}) where {T}
#     v[1] = t[1,1];
#     v[2] = t[2,2];
#     v[3] = t[3,3];
#     v[4] = t[1,2];
#     v[5] = t[2,1];
#     v[6] = t[2,3];
#     v[7] = t[3,2];
#     v[8] = t[1,3];
#     v[9] = t[3,1];
#     return v
# end

# """
#     strain9vto3x3t!(v::FVec{T}, t::FMat{T}) where {T}

# Convert a matrix of 3x3 strain components to a 9-component vector.

# The  strain components are in the order
#      ex, ey, ez, gxy, gyx, gyz, gzy, gxz, gzx
# """
# function strain9vto3x3t!(t::FMat{T}, v::FVec{T}) where {T}
#     t[1,1] = v[1];
#     t[2,2] = v[2];
#     t[3,3] = v[3];
#     t[1,2] = v[4];
#     t[2,1] = v[5];
#     t[2,3] = v[6];
#     t[3,2] = v[7];
#     t[1,3] = v[8];
#     t[3,1] = v[9];
#     return t
# end

# """
#     stress2x2to3v!(v::FVec{T}, t::FMat{T})

# Convert a symmetric matrix of 2x2 stress components to a 3-component vector.
# """
# function stress2x2to3v!(v::FVec{T}, t::FMat{T}) where {T}
#     v[1] = t[1,1];
#     v[2] = t[2,2];
#     v[3] = 0.5*(t[1,2] + t[2,1]);
#     return v
# end

# """
#     stress3vto2x2t!(t::FMat{T}, v::FVec{T})

# Convert a 3-vector to a  matrix of 2x2 stress components (symmetric tensor).
# """
# function stress3vto2x2t!(t::FMat{T}, v::FVec{T}) where {T}
#     t[1,1] = v[1];
#     t[2,2] = v[2];
#     t[1,2] = v[3];
#     t[2,1] = v[3];
#     return t
# end

# """
#     stress3vto3x3t!(t::FMat{T}, v::FVec{T}) where {T}

# Convert a 3-vector to a matrix of 3x3 stress components (symmetric tensor).
# """
# function stress3vto3x3t!(t::FMat{T}, v::FVec{T}) where {T}
#     t[1,1] = v[1];
#     t[2,2] = v[2];
#     t[1,2] = v[3];
#     t[2,1] = v[3];
#     return t
# end

# """
#     stress4vto3x3t!(t::FMat{T}, v::FVec{T}) where {T}

# Convert a 4-vector to a  matrix of 3x3 stress components (tensor).

# Convert a 4-vector to a *symmetric*
# matrix of 3x3 stress components (tensor).  This is
# conversion routine that would be useful for plane strain or
# axially symmetric conditions.
# The stress vector components need to be ordered as:
#     sigmax, sigmay, tauxy, sigmaz,
# which is the ordering used for the plane-strain model reduction.
# Therefore, for axially symmetric analysis the components need to be
# reordered, as from the constitutive equation they come out
# as sigmax, sigmay, sigmaz, tauxy.
# """
# function stress4vto3x3t!(t::FMat{T}, v::FVec{T}) where {T}
#     t[1,1] = v[1];
#     t[2,2] = v[2];
#     t[1,2] = v[3];
#     t[2,1] = v[3];
#     t[3,3] = v[4];
#     return t
# end

# """
#     stress6vto3x3t!(t::FMat{T}, v::FVec{T}) where {T}

# Convert a 6-vector to a  matrix of 3x3 stress components (symmetric tensor).
# """
# function stress6vto3x3t!(t::FMat{T}, v::FVec{T}) where {T}
#     t[1,1] = v[1];
#     t[2,2] = v[2];
#     t[3,3] = v[3];
#     t[1,2] = v[4];
#     t[2,1] = v[4];
#     t[1,3] = v[5];
#     t[3,1] = v[5];
#     t[3,2] = v[6];
#     t[2,3] = v[6];
#     return t
# end

# """
#     stress3x3tto6v!(v::FVec{T}, t::FMat{T}) where {T}

# Convert a matrix of 3x3 stress components to a 6-component vector.
# """
# function stress3x3tto6v!(v::FVec{T}, t::FMat{T}) where {T}
#     v[1] = t[1,1];
#     v[2] = t[2,2];
#     v[3] = t[3,3];
#     v[4] = 1/2.0*(t[1,2] + t[2,1]);
#     v[5] = 1/2.0*(t[1,3] + t[3,1]);
#     v[6] = 1/2.0*(t[3,2] + t[2,3]);
#     return v
# end

# """
#     strain9vto6v!(t::FVec{T}, v::FVec{T}) where {T}

# Convert a strain 9-vector to a  strain 6-vector components (tensor).
# """
# function strain9vto6v!(t::FVec{T}, v::FVec{T}) where {T}
#     t[1] = v[1];
#     t[2] = v[2];
#     t[3] = v[3];
#     t[4] = v[4]+v[5];
#     t[5] = v[8]+v[9];
#     t[6] = v[6]+v[7];
#     return t
# end

# """
#     strain6vto9v!(t::FVec{T}, v::FVec{T}) where {T}

# Convert a strain 6-vector to a  strain 9-vector components (tensor).

# The  strain components are in the order
#      ex, ey, ez, gxy/2, gxy/2, gyz/2, gyz/2, gxz/2, gxz/2
# """
# function strain6vto9v!(t::FVec{T}, v::FVec{T}) where {T}
#     t[1] = v[1];
#     t[2] = v[2];
#     t[3] = v[3];
#     t[4] = v[4]/2.;
#     t[5] = v[4]/2.;
#     t[6] = v[6]/2.;
#     t[7] = v[6]/2.;
#     t[8] = v[5]/2.;
#     t[9] = v[5]/2.;
#     return t
# end

# """
#     stress9vto6v!(t::FVec{T}, v::FVec{T}) where {T}

# Convert a stress 9-vector (tensor) to a  stress 6-vector components.

# The  stress components are in the order
#      sigx, sigy, sigz, tauxy, tauxy, tauyz, tauyz, tauxz, tauxz
# """
# function stress9vto6v!(t::FVec{T}, v::FVec{T}) where {T}
#     t[1] = v[1];
#     t[2] = v[2];
#     t[3] = v[3];
#     t[4] = v[4];
#     t[5] = v[8];
#     t[6] = v[6];
#     return t
# end

# """
#     stress6vto9v!(t::FVec{T}, v::FVec{T}) where {T}

# Convert a stress 6-vector to a  stress 9-vector components (tensor).

# The  stress components are in the order
#      sigx, sigy, sigz, tauxy, tauxy, tauyz, tauyz, tauxz, tauxz
# """
# function stress6vto9v!(t::FVec{T}, v::FVec{T}) where {T}
#     t[1] = v[1];
#     t[2] = v[2];
#     t[3] = v[3];
#     t[4] = v[4];
#     t[5] = v[4];
#     t[6] = v[6];
#     t[7] = v[6];
#     t[8] = v[5];
#     t[9] = v[5];
#     return t
# end

# function stressvectorrotation{MR<:DeforModelRed2DStress}(::Type{MR},
#   Rm::FMat{T})
#
#     # Rm = columns are components of 'bar' basis vectors on the 'plain'
#     #      basis vectors
#     #
#     # Calculate the rotation of the 'plain' stress vector to the
#     # 'bar' coordinate system given by the columns of the rotation matrix Rm.
#     #
#     # Example:
#     # The stress vector "stress" is given in the material coordinate
#     # system defined by the orientation matrix Rm. The following two
#     # transformations are equivalent:
#     #
#     #         t = stress_6v_to_3x3t (mat,stress);
#     #         t = (Rm*t*Rm');# in global coordinate system
#     #         t = (outputRm'*t*outputRm);# in output coordinate system
#     #         stress =stress_3x3t_to_6v (mat,t);# in output coordinate system
#     #
#     #        stress =mat.stress_vector_rotation(outputRm)...
#     #                  *mat.stress_vector_rotation(Rm')...
#     #                      *stress;# in output coordinate system
#
#     a11=Rm[1,1]; a12=Rm[1,2];
#     a21=Rm[2,1]; a22=Rm[2,2];
#     # switch self.reduction
#     #     case {'axisymm','strain'}
#     #         T =[[   a11^2,   a21^2, 0,         2*a11*a21]
#     #             [   a12^2,   a22^2, 0,         2*a12*a22]
#     #             [       0,       0, 1,                 0]
#     #             [ a11*a12, a21*a22, 0, a11*a22 + a12*a21]];
#     #             case 'stress'
#     T =[[   a11^2   a21^2       2*a11*a21]
#         [   a12^2   a22^2       2*a12*a22]
#         [ a11*a12 a21*a22 a11*a22+a12*a21]];
#     return T
# end
#
# function strainvectorrotation{MR<:DeforModelRed2DStress}(::Type{MR},
#   Rm::FMat{T})
#     # Calculate the rotation matrix for a strain vector.
#     #
#     #   Tbar = strain_vector_rotation(self,Rm)
#     #
#     # Rm = columns are components of 'bar' basis vectors on the 'plain'
#     #      basis vectors
#     #
#     # Calculate the rotation of the 'plain' strain vector to the
#     # 'bar' coordinate system given by the columns of the rotation matrix Rm.
#     #
#
#     a11=Rm[1,1]; a12=Rm[1,2];
#     a21=Rm[2,1]; a22=Rm[2,2];
#     # switch self.reduction
#     #     case {'axisymm','strain'}
#     #         Tbar =[ [     a11^2,     a21^2, 0,           a11*a21]
#     #                 [     a12^2,     a22^2, 0,           a12*a22]
#     #                 [         0,         0, 1,                 0]
#     #                 [ 2*a11*a12, 2*a21*a22, 0, a11*a22 + a12*a21]];
#     #   case 'stress'
#     Tbar =[ [     a11^2,     a21^2,           a11*a21]
#            [     a12^2,     a22^2,           a12*a22]
#            [ 2*a11*a12, 2*a21*a22, a11*a22 + a12*a21]];
#     return Tbar
# end
#
#
# ################################################################################
# # 2-D plane axially symmetric model
#
# function stressvectorrotation(::Type{MR},
#   Rm::FMat{T}) where {MR<:DeforModelRed2DAxisymm}
#
#     # Rm = columns are components of 'bar' basis vectors on the 'plain'
#     #      basis vectors
#     #
#     # Calculate the rotation of the 'plain' stress vector to the
#     # 'bar' coordinate system given by the columns of the rotation matrix Rm.
#     #
#     # Example:
#     # The stress vector "stress" is given in the material coordinate
#     # system defined by the orientation matrix Rm. The following two
#     # transformations are equivalent:
#     #
#     #         t = stress_6v_to_3x3t (mat,stress);
#     #         t = (Rm*t*Rm');# in global coordinate system
#     #         t = (outputRm'*t*outputRm);# in output coordinate system
#     #         stress =stress_3x3t_to_6v (mat,t);# in output coordinate system
#     #
#     #        stress =mat.stress_vector_rotation(outputRm)...
#     #                  *mat.stress_vector_rotation(Rm')...
#     #                      *stress;# in output coordinate system
#
#     a11=Rm[1,1]; a12=Rm[1,2];
#     a21=Rm[2,1]; a22=Rm[2,2];
#     T =[[   a11^2    a21^2  0          2*a11*a21]
#         [   a12^2    a22^2  0          2*a12*a22]
#         [       0        0  1                  0]
#         [ a11*a12  a21*a22  0  a11*a22 + a12*a21]];
#     return T
# end
#
# function strainvectorrotation(::Type{MR},
#   Rm::FMat{T}) where {MR<:DeforModelRed2DAxisymm}
#     # Calculate the rotation matrix for a strain vector.
#     #
#     #   Tbar = strain_vector_rotation(self,Rm)
#     #
#     # Rm = columns are components of 'bar' basis vectors on the 'plain'
#     #      basis vectors
#     #
#     # Calculate the rotation of the 'plain' strain vector to the
#     # 'bar' coordinate system given by the columns of the rotation matrix Rm.
#     #
#
#     a11=Rm[1,1]; a12=Rm[1,2];
#     a21=Rm[2,1]; a22=Rm[2,2];
#     Tbar =[ [     a11^2      a21^2  0            a11*a21]
#            [     a12^2      a22^2  0            a12*a22]
#            [         0          0  1                  0]
#            [ 2*a11*a12  2*a21*a22  0  a11*a22 + a12*a21]];
#     return Tbar
# end
#
# ################################################################################
# # Generic versions of rotations of stiffness and compliance matrices
#
# function rotatestiffness!(::Type{MR}, D::FMat{T},
#                                                            Rm::FMat{T}) where {MR<:DeforModelRed}
#     # Rotate constitutive stiffness matrix of the material.
#     #
#     #         function D=transform_stiffness(self,D,Rm)
#     #
#     # Rotate constitutive stiffness matrix of the material to the
#     # coordinate system given by the columns of the rotation matrix Rm.
#     T =stressvectorrotation(MR,Rm);
#     D = T*D*T';
#     return D
# end

#
# function rotatecompliance!(::Type{MR}, C::FMat{T},
#   Rm::FMat{T}) where {MR<:DeforModelRed}
#     # Rotate constitutive compliance matrix of the material.
#     #
#     #   C = rotate_compliance(self,C,Rm)
#     #
#     # Rotate constitutive compliance matrix of the material to the
#     # coordinate system given by the columns of the rotation matrix Rm.
#     Tbar =strainvectorrotation(MR,Rm);
#     C = Tbar*C*Tbar';
#     return C
# end
# function stressvectorrotation(::Type{DeforModelRed2DStrain}, Rm::FMat{T})
#
#     # Rm = columns are components of 'bar' basis vectors on the 'plain'
#     #      basis vectors
#     #
#     # Calculate the rotation of the 'plain' stress vector to the
#     # 'bar' coordinate system given by the columns of the rotation matrix Rm.
#     #
#     # Example:
#     # The stress vector "stress" is given in the material coordinate
#     # system defined by the orientation matrix Rm. The following two
#     # transformations are equivalent:
#     #
#     #         t = stress_6v_to_3x3t (mat,stress);
#     #         t = (Rm*t*Rm');# in global coordinate system
#     #         t = (outputRm'*t*outputRm);# in output coordinate system
#     #         stress =stress_3x3t_to_6v (mat,t);# in output coordinate system
#     #
#     #        stress =mat.stress_vector_rotation(outputRm)...
#     #                  *mat.stress_vector_rotation(Rm')...
#     #                      *stress;# in output coordinate system
#
#     a11=Rm[1,1]; a12=Rm[1,2];
#     a21=Rm[2,1]; a22=Rm[2,2];
#     # sigmax, sigmay, sigmaz, tauxy
#     T =[[   a11^2    a21^2  0          2*a11*a21]
#         [   a12^2    a22^2  0          2*a12*a22]
#         [       0        0  1                  0]
#         [ a11*a12  a21*a22  0  a11*a22 + a12*a21]];
#     return T
# end
#
# function strainvectorrotation(::Type{DeforModelRed2DStrain},  Rm::FMat{T})
#     # Calculate the rotation matrix for a strain vector.
#     #
#     #   Tbar = strain_vector_rotation(self,Rm)
#     #
#     # Rm = columns are components of 'bar' basis vectors on the 'plain'
#     #      basis vectors
#     #
#     # Calculate the rotation of the 'plain' strain vector to the
#     # 'bar' coordinate system given by the columns of the rotation matrix Rm.
#     #
#
#     a11=Rm[1,1]; a12=Rm[1,2];
#     a21=Rm[2,1]; a22=Rm[2,2];
#     Tbar =[ [     a11^2      a21^2  0            a11*a21]
#            [     a12^2      a22^2  0            a12*a22]
#            [         0          0  1                  0]
#            [ 2*a11*a12  2*a21*a22  0  a11*a22 + a12*a21]];
#     return Tbar
# end
# function strainvectorrotation(::Type{DeforModelRed3D}, Rm::FMat{T})
#     # Calculate the rotation matrix for a strain vector.
#     #
#     #   Tbar = strain_vector_rotation(self,Rm)
#     #
#     # Rm = columns are components of 'bar' basis vectors on the 'plain'
#     #      basis vectors
#     #
#     # Calculate the rotation of the 'plain' strain vector to the
#     # 'bar' coordinate system given by the columns of the rotation matrix Rm.
#     #
#
#     a11=Rm[1,1]; a12=Rm[1,2]; a13=Rm[1,3];
#     a21=Rm[2,1]; a22=Rm[2,2]; a23=Rm[2,3];
#     a31=Rm[3,1]; a32=Rm[3,2]; a33=Rm[3,3];
#     Tbar =[
#            [     a11^2     a21^2     a31^2           a11*a21           a11*a31           a21*a31]
#            [     a12^2     a22^2     a32^2           a12*a22           a12*a32           a22*a32]
#            [     a13^2     a23^2     a33^2           a13*a23           a13*a33           a23*a33]
#            [ 2*a11*a12 2*a21*a22 2*a31*a32 a11*a22 + a12*a21 a11*a32 + a12*a31 a21*a32 + a22*a31]
#            [ 2*a11*a13 2*a21*a23 2*a31*a33 a11*a23 + a13*a21 a11*a33 + a13*a31 a21*a33 + a23*a31]
#            [ 2*a12*a13 2*a22*a23 2*a32*a33 a12*a23 + a13*a22 a12*a33 + a13*a32 a22*a33 + a23*a32]];
#     return Tbar
# end


text:
export onehot, thermometer

include("onehot.jl")
include("thermometer.jl")


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm

V,EV = Lar.randomcuboids(250, .2)
V = Plasm.normalize(V,flag=true)
Plasm.view(Plasm.numbering(.05)((V,[[[k] for k=1:size(V,2)], EV])))

W = convert(Lar.Points, V')
cop_EV = Lar.coboundary_0(EV::Lar.Cells)
cop_EW = convert(Lar.ChainOp, cop_EV)
V, copEV, copFE = Lar.Arrangement.planar_arrangement(W::Lar.Points, cop_EW::Lar.ChainOp)

triangulated_faces = Lar.triangulate2D(V, [copEV, copFE])
FVs = convert(Array{Lar.Cells}, triangulated_faces)
V = convert(Lar.Points, V')
Plasm.viewcolor(V::Lar.Points, FVs::Array{Lar.Cells})

W, copEV, copFE = Lar.Arrangement.planar_arrangement(W::Lar.Points, cop_EW::Lar.ChainOp)
EVs = Lar.FV2EVs(copEV, copFE) # polygonal face fragments
V = convert(Lar.Points, W')
Plasm.viewcolor(V::Lar.Points, EVs::Array{Lar.Cells})

model = V,EVs
Plasm.view(Plasm.lar_exploded(model)(1.2,1.2,1.2))


text:
module GeometricFlux

using Statistics: mean
using SparseArrays: SparseMatrixCSC
using LinearAlgebra: Adjoint, norm

using FillArrays: Fill
using Flux
using Flux: glorot_uniform, leakyrelu, GRUCell
using Flux: @functor
using GraphSignals
using LightGraphs
using Requires
using ScatterNNlib
using Zygote
using ZygoteRules

import Flux: maxpool, meanpool

export
    # layers/gn
    GraphNet,
    update_edge,
    update_vertex,
    update_global,
    update_batch_edge,
    update_batch_vertex,
    aggregate_neighbors,
    aggregate_edges,
    aggregate_vertices,
    propagate,

    # layers/msgpass
    MessagePassing,
    message,
    update,

    # layers/conv
    GCNConv,
    ChebConv,
    GraphConv,
    GATConv,
    GatedGraphConv,
    EdgeConv,
    message,
    update,
    propagate,

    # layer/pool
    GlobalPool,
    LocalPool,
    TopKPool,

    # models
    GAE,
    VGAE,
    InnerProductDecoder,
    VariationalEncoder,
    summarize,
    sample,

    # layer/selector
    FeatureSelector,

    # operations/pool
    sumpool,
    subpool,
    prodpool,
    divpool,
    maxpool,
    minpool,
    meanpool,
    pool,

    # graph/index
    generate_cluster,

    # utils
    topk_index

const IntOrTuple = Union{Integer,Tuple}

include("pool.jl")

include("graph/index.jl")

include("utils.jl")

include("layers/gn.jl")
include("layers/msgpass.jl")

include("layers/conv.jl")
include("layers/pool.jl")
include("models.jl")
include("layers/selector.jl")

include("graph/simplegraphs.jl")

function __init__()
    @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" begin
        using CUDA
        include("cuda/msgpass.jl")
        include("cuda/conv.jl")
        include("cuda/pool.jl")
    end
    @require SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622" begin
        include("graph/weightedgraphs.jl")
    end
    @require MetaGraphs = "626554b9-1ddb-594c-aa3c-2596fe9399a5" begin
        include("graph/metagraphs.jl")
    end
end

end


text:
@testset "1017.convert-to-base-2.jl" begin
    @test base_neg2(3) == "111"
    @test base_neg2(2) == "110"
    @test base_neg2(4) == "100"
    @test base_neg2(12524) == "111000100111100"
end


text:
# Define the solver type
abstract type Solver end

function integrate(a!::Function, b!::Function, x0::AbstractVector{T}, t::AbstractVector{<:Real}, solver::Solver, m::Int, seed::Integer;
                params...) where {T<:Number}
    noise = Noise(m;seed=seed)
    out = [[], []] #Instantiate for saving the x and W respectively
    callback = default_callback(out)
    integrate(a!, b!, x0, t, noise, solver, callback; params...)
    return out
end

max_stars_repo_name: onoderat/SDEs

max_issues_repo_name: onoderat/SDEs

max_forks_repo_name: onoderat/SDEs


text:
"""
Julia interface to a local verison of the empca
(https://github.com/sbailey/empca) python package
"""
__precompile__()
module EMPCA
    using PyCall
    const empca = PyNULL()
    function __init__()
        pyimport_conda("scipy", "scipy")
        pushfirst!(PyVector(pyimport("sys")."path"), @__DIR__)
        copy!(empca, pyimport("empca"))
        @assert empca != PyNULL()
    end
    export empca
end


text:
using Flux, Flux.Tracker, CuArrays, Base.Test
using Flux: gpu

info("Testing Flux/GPU")

@testset "CuArrays" begin

CuArrays.allowscalar(false)

x = param(randn(5, 5))
cx = gpu(x)
@test cx isa TrackedArray && cx.data isa CuArray

x = Flux.onehotbatch([1, 2, 3], 1:3)
cx = gpu(x)
@test cx isa Flux.OneHotMatrix && cx.data isa CuArray

m = Chain(Dense(10, 5, tanh), Dense(5, 2), softmax)
cm = gpu(m)

@test all(p isa TrackedArray && p.data isa CuArray for p in params(cm))
@test cm(gpu(rand(10, 10))) isa TrackedArray{Float32,2,CuArray{Float32,2}}

x = [1,2,3]
cx = gpu(x)
@test Flux.crossentropy(x,x) ≈ Flux.crossentropy(cx,cx)

# Fails in Pkg.test ffs
# c = gpu(Conv((2,2),3=>4))
# l = c(gpu(rand(10,10,3,2)))
# Flux.back!(sum(l))

end

CuArrays.cudnn_available() && include("cudnn.jl")


text:
function DrawTri!(T)
        plot!([T[1,1],T[2,1]],[T[1,2],T[2,2]],label="")
        plot!([T[2,1],T[3,1]],[T[2,2],T[3,2]],label="")
        plot!([T[3,1],T[1,1]],[T[3,2],T[1,2]],label="")
end

function Ellipse(n;a=1,b=1,max=2π,θ=0,center=[0.0,0.0])
    x = [a*cos(t)*cos(θ) - b*sin(t)*sin(θ) for t in 0:(2π/n):2π] .+ center[1]
    y = [a*cos(t)*sin(θ) + b*sin(t)*cos(θ) for t in 0:(2π/n):2π] .+ center[2]
    return cat(x,y,dims=2)
end

function PackingFraction(X,ParticleArea)
    return (size(X,2)*ParticleArea) ./ Alphashapes.AlphaShapeVolume(AlphaShape(X,FindAlpha(X)))
end

max_issues_repo_name: harveydevereux/CUDA-Whirligigs

max_forks_repo_name: harveydevereux/CUDA-Whirligigs


text:
export
    LineSegment,
    parallel

immutable LineSegment
    A::VecE2
    B::VecE2
end
Base.:+(seg::LineSegment, V::VecE2) = LineSegment(seg.A + V, seg.B + V)
Base.:-(seg::LineSegment, V::VecE2) = LineSegment(seg.A - V, seg.B - V)
Base.convert(::Type{Line}, seg::LineSegment) = Line(seg.A, seg.B)
Base.convert(::Type{LineSegment}, line::Line) = LineSegment(line.A, line.B)

get_polar_angle(seg::LineSegment) = mod2pi(atan2(seg.B.y - seg.A.y, seg.B.x - seg.A.x))

"""
The distance between the line segment and the point P
"""
function get_distance(seg::LineSegment, P::VecE2)

    ab = seg.B - seg.A
    pb = P - seg.A

    denom = abs2(ab)
    if denom == 0.0
        return 0.0
    end

    r = dot(ab, pb)/denom

    if r ≤ 0.0
        abs(P - seg.A)
    elseif r ≥ 1.0
        abs(P - seg.B)
    else
        abs(P - (seg.A + r*ab))
    end
end

"""
What side of the line you are on, based on A → B
"""
get_side(seg::LineSegment, p::VecE2) = sign((seg.B.x-seg.A.x) * (p.y-seg.A.y) - (seg.B.y-seg.A.y) * (p.x-seg.A.x))

"""
The angular distance between the two line segments
"""
function angledist(segA::LineSegment, segB::LineSegment)
    u = segA.B - segA.A
    v = segB.B - segB.A
    sqdenom = abs2(u)*abs2(v)
    if isapprox(sqdenom, 0.0, atol=1e-10)
        return NaN
    end
    return acos(dot(u,v) / sqrt(sqdenom))
end

"""
True if the two segments are parallel
"""
function parallel(segA::LineSegment, segB::LineSegment, ε::Float64=1e-10)
    θ = angledist(segA, segB)
    return isapprox(θ, 0.0, atol=ε)
end


text:
export lorenz96!, generate_lorenz96, spin_lorenz96, setup_lorenz96, benchmark_lorenz96, setup_order, benchmark_srmf_lorenz96, benchmark_sadaptivermf_lorenz96

"""
    lorenz96!(du,u,p,t)

Compute in-place the right-hand-side of the Lorenz-96 system for a state `u` at time `t`,
and store it in `du`. `p` is vector of user-defined parameters.
"""
function lorenz96!(du,u,p,t)
    F = 8.0
    n = size(u,1)
    du[1] = (u[2]-u[end-1])*u[end] - u[1] + F
    du[2] = (u[3]-u[end])*u[1] - u[2] + F
    du[end] = (u[1] - u[end-2])*u[end-1] - u[end] + F

    @inbounds for i=3:n-1
        du[i] = (u[i+1] - u[i-2])*u[i-1] - u[i] + F
    end
    return du
end

function generate_lorenz96(model::Model, x0::Array{Float64,1}, J::Int64)

    @assert model.Nx == size(x0,1) "Error dimension of the input"
    xt = zeros(model.Nx,J)
    x = deepcopy(x0)
    yt = zeros(model.Ny,J)
    tt = zeros(J)

    t0 = 0.0

    step = ceil(Int, model.Δtobs/model.Δtdyn)
    tspan = (t0, t0 + model.Δtobs)
    prob = ODEProblem(model.F.f,x,tspan)

    for i=1:J
    	# Run dynamics and save results
    	tspan = (t0 + (i-1)*model.Δtobs, t0 + i*model.Δtobs)
    	prob = remake(prob, tspan = tspan)

    	sol = solve(prob, Tsit5(), dt = model.Δtdyn, adaptive = false,
		            dense = false, save_everystep = false)
    	x .= deepcopy(sol.u[end])
    	# for j=1:step
    	# 	t = t0 + (i-1)*algo.Δtobs+(j-1)*algo.Δtdyn
        # 	_, x = model.f(t+(i-1)*model.Δtdyn, x)
    	# end
    	model.ϵx(x)

    	# Collect observations
    	tt[i] = deepcopy(i*model.Δtobs)
    	xt[:,i] = deepcopy(x)
		yt[:,i] = deepcopy(model.F.h(x, tt[i]))
		yt[:,i] .+= model.ϵy.m + model.ϵy.σ*randn(model.Ny)
    end
    	return SyntheticData(tt, model.Δtdyn, x0, xt, yt)
end

function spin_lorenz96(model::Model, data::SyntheticData, Ne::Int64, path::String)
	@assert path[1]=='/' && path[end]=='/' "This function expects a / at the extremities of path"

	# Set initial condition
	X = zeros(model.Ny + model.Nx, Ne)
	X[model.Ny+1:model.Ny+model.Nx,:] .= rand(model.π0, Ne)#sqrt(model.C0)*randn(model.Nx, Ne) .+ model.m0

	J = model.Tspinup
	t0 = 0.0
	F = model.F
	enkf = StochEnKF(x->x, model.ϵy, model.Δtdyn, model.Δtobs, false, false)

	statehist = seqassim(F, data, J, model.ϵx, enkf, X, model.Ny, model.Nx, t0)

    save(path*"set_up_Ne"*string(Ne)*".jld", "state", statehist, "Ne", Ne, "x0", data.x0, "tt", data.tt, "xt", data.xt, "yt", data.yt)

    # return statehist
	_,_,rmse_mean,_ = metric_hist(rmse, data.xt[:,1:J], statehist[2:end])
	println("Ne "*string(Ne)* " RMSE: "*string(rmse_mean))
	# Save data
	save(path*"set_up_Ne"*string(Ne)*".jld", "X", statehist[end], "Ne", Ne, "x0", data.x0, "tt", data.tt, "xt", data.xt, "yt", data.yt)
end

function setup_lorenz96(path::String, Ne_array::Array{Int64,1})
    Nx = 40
    Ny = 20

    # The state is measured every fourth
    Δ = 2
    yidx = 1:Δ:Nx

    Δtdyn = 0.01
    Δtobs = 0.4

    σx = 1e-6#1e-2#1e-6#1e-2
    σy = sqrt(0.5)#1e-6#2.0

    ϵx = AdditiveInflation(Nx, zeros(Nx), σx)
    ϵy = AdditiveInflation(Ny, zeros(Ny), σy)
    tspinup = 200.0
    Tspinup = 2000
    tmetric = 400.0
    Tmetric = 4000
    t0 = 0.0
    tf = 600.0
    Tf = 6000

    Tburn = 2000
    Tstep = 4000
    Tspinup = 2000

    f = lorenz96!
    h(x, t) = x[yidx]
    # Create a local version of the observation operator
    h(x, idx, t) = x[idx]
	F = StateSpace(f, h)

    model = Model(Nx, Ny, Δtdyn, Δtobs, ϵx, ϵy, MvNormal(zeros(Nx), Matrix(1.0*I, Nx, Nx)), Tburn, Tstep, Tspinup, F);

    # Set initial condition
    x0 = rand(model.π0)
    # x0 = [0.645811242103507;  -1.465126216973632;   0.385227725149822];

    # Run dynamics and generate data
    data = generate_lorenz96(model, x0, model.Tspinup+model.Tstep);

    for Ne in Ne_array
        spin_lorenz96(model, data, Ne, path)
    end

    return model, data
end

function benchmark_lorenz96(model::Model, data::SyntheticData, path::String, Ne_array::Array{Int64,1},
	                        β_array::Array{Float64,1}, Lrad_array::Array{Float64,1})
@assert path[1]=='/' && path[end]=='/' "This function expects a / at the extremities of path"

#Store all the metric per number of ensemble members
metric_list = []

@showprogress for Ne in Ne_array
    metric_Ne = Metrics[]
    for β in β_array
		for Lrad in Lrad_array
		    @show Ne, β, Lrad
		    # Load file
		    X0 = load(path*"set_up_Ne"*string(Ne)*".jld", "X")

		    X = zeros(model.Ny + model.Nx, Ne)
		    X[model.Ny+1:model.Ny+model.Nx,:] .= copy(X0)
		    J = model.Tstep
		    t0 = model.Tspinup*model.Δtobs
		    F = model.F

			Δ = 2
			yidx = 1:Δ:model.Nx
			idx = vcat(collect(1:length(yidx))', collect(yidx)')

			enkf = SeqStochEnKF(x->x, F.h, MultiplicativeInflation(β), model.ϵy, model.Ny, model.Nx,
			                    model.Δtdyn, model.Δtobs,
			                    idx, zeros(model.Nx+1, Ne), false, true)
		    # Use this multi-additive inflation
		    ϵx = MultiAddInflation(model.Nx, β, zeros(model.Nx), model.ϵx.Σ, model.ϵx.σ)

			# Create Localization structure
		    Gxy(i,j) = periodicmetric!(i,yidx[j], model.Nx)
		    Gyy(i,j) = periodicmetric!(yidx[i],yidx[j], model.Nx)
		    Loc = Localization(Lrad, Gxy, Gyy)

			@time statehist = seqassim(F, data, J, model.ϵx, enkf, X, model.Ny, model.Nx, t0, Loc);

		    metric = post_process(data, model, J, statehist)
		    push!(metric_Ne, deepcopy(metric))
			println("Ne "*string(Ne)*"& β "*string(β)*"& Lrad "*string(Lrad)*" RMSE: "*string(metric.rmse_mean))
		end
	end
	push!(metric_list, deepcopy(metric_Ne))
end

return metric_list
end

# Nx size of the state
# off_p contains the order for the off-diagonal component
# off_rad contains the radius for the localization of the map
# nonid_rad contains the number of non id components
# dist matrix of the distances between the variables
function setup_order(Nx, diagobs_p, diagunobs_p, off_p, off_rad, nonid_rad, dist::Array{Float64,2})

	dist = convert(Array{Int64,2}, dist)
    perm = sortperm(view(dist,:,1))

    dist_to_order = fill(-1, Nx)
    # fill!(view(dist_to_order,1:ceil(Int64,off_rad)), off_p)
    fill!(view(dist_to_order,1:ceil(Int64,min(off_rad, Nx))), off_p)
    order = [fill(-1,i) for i=1:Nx+1]
    nonid = Int64[]
    for i=1:Nx
        node_i = perm[i]
        dist_i1 = dist[1, node_i]

        # Vector to store order, offset of 1
        # because we have one scalar observation
        order_i = zeros(Int64, i-1)#view(order,i+1)
        if dist_i1 <= nonid_rad
            if i>1
            for j=1:i-1

                # extract distance to node i
                node_j = perm[j]
                dist_ij = dist[node_i, node_j]

                # compute order for variable based on distance
                order_i[j] = dist_to_order[dist_ij]
            end
            end
            # Define order for diagonal TM components and observation component
            if i==1 # observed mode
                order[2][2] = diagobs_p
                order[2][1] = off_p
            else # unobserved mode
                order[i+1][2:i] .= deepcopy(order_i)
                order[i+1][i+1] = diagunobs_p
                order[i+1][1] = -1
            end
        end
    end
    return order
end

# Routine to benchmark the stochastic radial map filter
function benchmark_srmf_lorenz96(model::Model, data::SyntheticData, path::String, Ne_array::Array{Int64,1}, β_array::Array{Float64,1},
	                             Lrad_array::Array{Float64,1}, nonid_array::Array{Float64,1}, p::Int64)
@assert path[1]=='/' && path[end]=='/' "This function expects a / at the extremities of path"

#Store all the metric per number of ensemble members
metric_hist = []

@showprogress for Ne in Ne_array
    metric_Ne = Metrics[]
    for β in β_array
		for Lrad in Lrad_array
			for nonid_rad in nonid_array
					@show Ne, β, Lrad, nonid_rad
				    # Load file
				    X0 = load(path*"set_up_Ne"*string(Ne)*".jld", "X")

				    X = zeros(model.Ny + model.Nx, Ne)
				    X[model.Ny+1:model.Ny+model.Nx,:] .= copy(X0)

				    J = model.Tstep
				    t0 = model.Tspinup*model.Δtobs
				    F = model.F

					# parameters of the radial map
					γ = 2.0
					λ = 0.01
					δ = 1e-8
					κ = 4.0

					dist = Float64.(metric_lorenz(model.Nx))
					idx = vcat(collect(1:model.Ny)',collect(1:2:model.Nx)')
					order = setup_order(model.Nx, p, 0, p, Lrad, nonid_rad, dist)

					# Define index of the measurement

					smf = SparseRadialSMF(x->x, F.h, β, model.ϵy, order, γ, λ, δ, κ,
				                      model.Ny, model.Nx, Ne,
				                      model.Δtdyn, model.Δtobs,
				                      dist, idx; islocalized = true)

				    @time statehist = seqassim(F, data, J, model.ϵx, smf, X, model.Ny, model.Nx, t0);

				    metric = post_process(data, model, J, statehist)
				    push!(metric_Ne, deepcopy(metric))
				    println("Ne "*string(Ne)*" & β "*string(β)*" & Lrad "*string(Lrad)*" & nonid_rad "*string(nonid_rad)*" & RMSE: "*string(metric.rmse_mean))
				end
			end
	    end
	    push!(metric_hist, deepcopy(metric_Ne))
	end
	return metric_hist
end

# Routine to benchmark the stochastic adaptive radial map filter
function benchmark_sadaptivermf_lorenz96(model::Model, data::SyntheticData, path::String, Ne_array::Array{Int64,1},
	                                     β_array::Array{Float64,1}, nonid_array::Array{Float64,1}, p::Int64)
@assert path[1]=='/' && path[end]=='/' "This function expects a / at the extremities of path"

# Perform a free-run of the model to identify the structure

xfreerun = rand(model.π0)
tfreerun = 10000
Tfreerun = ceil(Int64, tfreerun/model.Δtobs)
data_freerun = generate_lorenz96(model, xfreerun, Tfreerun)

Tdiscard = 1000
Δdiscard = 10

#Store all the metric per number of ensemble members
metric_hist = []

@showprogress for Ne in Ne_array
    metric_Ne = Metrics[]
    for β in β_array
		for nonid_rad in nonid_array

			γ = 2.0
			λ = 0.0
			δ = 1e-8
			κ = 10.0

			J = model.Tstep
			t0 = model.Tspinup*model.Δtobs
			F = model.F
			dist = Float64.(metric_lorenz(model.Nx))
			idx = vcat(collect(1:model.Ny)',collect(1:model.Ny)')

			Slist = SparseRadialMap[]
				for i=1:model.Ny
					idx1, idx2 = idx[:,i]
					perm = sortperm(view(dist,:,idx2))
					Xi = vcat(data_freerun.yt[i:i,Tdiscard:Δdiscard:end], data_freerun.xt[perm, Tdiscard:Δdiscard:end])
					Si = SparseRadialMap(Xi, -1; λ = λ, δ = δ, γ = γ)
					poff = p
					order = [p; 0; 0]
					# maxfeatures = ceil(Int64, (sqrt(Ne)-(findmin(order)[2]+1))/(poff+1))

					optimize(Si, Xi, poff, order, 2; apply_rescaling = true, start = 2)#, maxfeatures = maxfeatures)
					push!(Slist, Si)
				end
		    # Load file
		    X0 = load(path*"set_up_Ne"*string(Ne)*".jld", "X")

		    X = zeros(model.Ny + model.Nx, Ne)
		    X[model.Ny+1:model.Ny+model.Nx,:] .= copy(X0)

		    smf = AdaptiveSparseRadialSMF(x->x, F.h, β, model.ϵy, Slist,
		                        model.Ny, model.Nx,
		                        model.Δtdyn, model.Δtobs, Inf,
		                        dist, idx, zeros(model.Nx+1, Ne), false, true)
		    @time statehist = seqassim(F, data, J, model.ϵx, smf, X, model.Ny, model.Nx, t0);

		    metric = post_process(data, model, J, statehist)
		    push!(metric_Ne, deepcopy(metric))
		    println("Ne "*string(Ne)*"& β "*string(β)*" & nonid_rad "*string(nonid_rad)*" RMSE: "*string(metric.rmse_mean))
		    end
		end
		push!(metric_hist, deepcopy(metric_Ne))
	end

return metric_hist
end


text:
@testset "Online" begin

CI_USER_CFG = aws_config()
# do not set this variable in CI; it should be versioned with the code
# this is for locally overriding the stack used in testing
TEST_STACK_NAME = get(ENV, "CLOUDWATCHLOGSJL_STACK_NAME", "CloudWatchLogs-jl-00015")
TEST_RESOURCE_PREFIX = "pubci-$TEST_STACK_NAME-cwl-test"
TEST_LOG_GROUP = "$TEST_RESOURCE_PREFIX-group"
FORBIDDEN_LOG_GROUP = "$TEST_RESOURCE_PREFIX-group-forbidden"
FORBIDDEN_GROUP_LOG_STREAM = "$TEST_RESOURCE_PREFIX-group-forbidden-stream"
BAD_STREAM_LOG_GROUP = "$TEST_RESOURCE_PREFIX-group-badstream"
FORBIDDEN_LOG_STREAM = "$TEST_RESOURCE_PREFIX-stream-forbidden"
TEST_ROLE = stack_output(CI_USER_CFG, TEST_STACK_NAME)["LogTestRoleArn"]
CFG = aws_config(creds=assume_role(CI_USER_CFG, TEST_ROLE; DurationSeconds=7200))
LOG_RUN_ID = uuid1()

new_stream = let
    counter = 1

    function new_stream(category::AbstractString)
        stream_name = @sprintf "pubci-%s-%03d-%s" category counter LOG_RUN_ID
        counter += 1
        return stream_name
    end
end

new_group = let
    counter = 1

    function new_group(category::AbstractString)
        stream_name = @sprintf "pubci-%s-%03d-%s" category counter LOG_RUN_ID
        counter += 1
        return stream_name
    end
end

@testset "Create/delete groups and streams" begin
    @testset "Named group" begin
        group_name = new_group("create_group")
        @test create_group(CFG, group_name; tags=Dict("Temporary"=>"true")) == group_name

        response = logs(
            CFG,
            "DescribeLogGroups";
            logGroupNamePrefix=group_name,
            limit=1,
        )

        groups = response["logGroups"]

        @test !isempty(groups)
        @test groups[1]["logGroupName"] == group_name

        delete_group(CFG, group_name)

        response = logs(
            CFG,
            "DescribeLogGroups";
            logGroupNamePrefix=group_name,
            limit=1,
        )

        groups = response["logGroups"]

        @test isempty(groups) || groups[1]["logGroupName"] != group_name
    end

    @testset "Named group no tags" begin
        group_name = new_group("create_group_no_tags")
        @test create_group(CFG, group_name) == group_name

        response = logs(
            CFG,
            "DescribeLogGroups";
            logGroupNamePrefix=group_name,
            limit=1,
        )

        groups = response["logGroups"]

        @test !isempty(groups)
        @test groups[1]["logGroupName"] == group_name

        delete_group(CFG, group_name)

        response = logs(
            CFG,
            "DescribeLogGroups";
            logGroupNamePrefix=group_name,
            limit=1,
        )

        groups = response["logGroups"]

        @test isempty(groups) || groups[1]["logGroupName"] != group_name
    end

    @testset "Unnamed group" begin
        group_name = create_group(CFG; tags=Dict("Temporary"=>"true"))

        response = logs(
            CFG,
            "DescribeLogGroups";
            logGroupNamePrefix=group_name,
            limit=1,
        )

        groups = response["logGroups"]

        @test !isempty(groups)
        @test groups[1]["logGroupName"] == group_name

        delete_group(CFG, group_name)

        response = logs(
            CFG,
            "DescribeLogGroups";
            logGroupNamePrefix=group_name,
            limit=1,
        )

        groups = response["logGroups"]

        @test isempty(groups) || groups[1]["logGroupName"] != group_name
    end

    @testset "Not allowed" begin
        @test_throws AWSException create_group(CFG, "delta∆")  # invalid characters
        @test_throws AWSException delete_group(CFG, TEST_LOG_GROUP)  # explicitly forbidden
    end

    @testset "Named stream" begin
        stream_name = new_stream("create_stream")
        @test create_stream(CFG, TEST_LOG_GROUP, stream_name) == stream_name

        response = logs(
            CFG,
            "DescribeLogStreams";
            logGroupName=TEST_LOG_GROUP,
            logStreamNamePrefix=stream_name,
            orderBy="LogStreamName",  # orderBy and limit will ensure we get just the one
            limit=1,                  # matching result
        )

        streams = response["logStreams"]

        @test !isempty(streams)
        @test streams[1]["logStreamName"] == stream_name

        delete_stream(CFG, TEST_LOG_GROUP, stream_name)

        response = logs(
            CFG,
            "DescribeLogStreams";
            logGroupName=TEST_LOG_GROUP,
            logStreamNamePrefix=stream_name,
            orderBy="LogStreamName",  # orderBy and limit will ensure we get just the one
            limit=1,                  # matching result
        )

        streams = response["logStreams"]

        @test isempty(streams) || streams[1]["logStreamName"] != stream_name
    end

    @testset "Unnamed stream" begin
        stream_name = create_stream(CFG, TEST_LOG_GROUP)

        response = logs(
            CFG,
            "DescribeLogStreams";
            logGroupName=TEST_LOG_GROUP,
            logStreamNamePrefix=stream_name,
            orderBy="LogStreamName",  # orderBy and limit will ensure we get just the one
            limit=1,                  # matching result
        )

        streams = response["logStreams"]

        @test !isempty(streams)
        @test streams[1]["logStreamName"] == stream_name

        delete_stream(CFG, TEST_LOG_GROUP, stream_name)

        response = logs(
            CFG,
            "DescribeLogStreams";
            logGroupName=TEST_LOG_GROUP,
            logStreamNamePrefix=stream_name,
            orderBy="LogStreamName",  # orderBy and limit will ensure we get just the one
            limit=1,                  # matching result
        )

        streams = response["logStreams"]

        @test isempty(streams) || streams[1]["logStreamName"] != stream_name
    end

    @testset "Not allowed" begin
        @test_throws AWSException create_stream(CFG, FORBIDDEN_LOG_GROUP)
        @test_throws AWSException delete_stream(CFG, FORBIDDEN_LOG_GROUP, FORBIDDEN_GROUP_LOG_STREAM)
    end
end

@testset "CloudWatchLogStream" begin
    @testset "Normal log submission" begin
        start_time = CloudWatchLogs.unix_timestamp_ms()
        stream_name = new_stream("stream_type")
        @test create_stream(CFG, TEST_LOG_GROUP, stream_name) == stream_name

        stream = CloudWatchLogStream(CFG, TEST_LOG_GROUP, stream_name)
        @test submit_log(stream, LogEvent("Hello AWS")) == 1
        @test submit_logs(stream, LogEvent.(["Second log", "Third log"])) == 2

        sleep(2)  # wait until AWS has injested the logs; this may or may not be enough
        response = logs(
            CFG,
            "GetLogEvents";
            logGroupName=TEST_LOG_GROUP,
            logStreamName=stream_name,
            startFromHead=true,
        )

        time_range = (start_time - 10):(CloudWatchLogs.unix_timestamp_ms() + 10)

        @test length(response["events"]) == 3
        messages = map(response["events"]) do event
            @test round(Int64, event["timestamp"]) in time_range
            event["message"]
        end

        @test messages == ["Hello AWS", "Second log", "Third log"]
        delete_stream(CFG, TEST_LOG_GROUP, stream_name)
    end

    @testset "Not allowed" begin
        @test_throws AWSException CloudWatchLogStream(CFG, FORBIDDEN_LOG_GROUP, FORBIDDEN_GROUP_LOG_STREAM)

        stream = CloudWatchLogStream(CFG, BAD_STREAM_LOG_GROUP, FORBIDDEN_LOG_STREAM, nothing)
        @test_throws AWSException submit_log(stream, LogEvent("Foo"))
    end

    @testset "Too many logs" begin
        start_time = CloudWatchLogs.unix_timestamp_ms()

        stream = CloudWatchLogStream(
            CFG,
            TEST_LOG_GROUP,
            create_stream(CFG, TEST_LOG_GROUP, new_stream("too_many_logs")),
        )

        events = map(Iterators.take(Iterators.countfrom(start_time), 10001)) do ts
            LogEvent("A", ts)
        end

        @test_throws LOGGER LogSubmissionException submit_logs(stream, events)
    end

    @testset "Logs too big" begin
        stream = CloudWatchLogStream(
            CFG,
            TEST_LOG_GROUP,
            create_stream(CFG, TEST_LOG_GROUP, new_stream("logs_too_big")),
        )

        event_size = CloudWatchLogs.MAX_EVENT_SIZE - 26
        events = map(1:(div(CloudWatchLogs.MAX_BATCH_SIZE, event_size) + 1)) do i
            LogEvent("A" ^ event_size)
        end
        @test_throws LOGGER LogSubmissionException submit_logs(stream, events)
    end

    @testset "Logs too spread" begin
        stream = CloudWatchLogStream(
            CFG,
            TEST_LOG_GROUP,
            create_stream(CFG, TEST_LOG_GROUP, new_stream("logs_too_spread")),
        )

        last_time = Dates.now(tz"UTC")
        first_time = last_time - Hour(25)

        events = [LogEvent("First", first_time), LogEvent("Last", last_time)]

        @test_throws LOGGER LogSubmissionException submit_logs(stream, events)
    end

    @testset "Invalid sequence token" begin
        stream = CloudWatchLogStream(
            CFG,
            TEST_LOG_GROUP,
            create_stream(CFG, TEST_LOG_GROUP, new_stream("invalid_token")),
        )

        @test submit_log(stream, LogEvent("Foo")) == 1
        CloudWatchLogs.update_sequence_token!(stream, "oops_invalid")
        setlevel!(LOGGER, "debug") do
            @test_log LOGGER "debug" "InvalidSequenceTokenException" begin
                submit_log(stream, LogEvent("Second time's the charm"))
            end
        end
    end

    @testset "Unsorted logs" begin
        stream_name = new_stream("unsorted")
        stream = CloudWatchLogStream(
            CFG,
            TEST_LOG_GROUP,
            create_stream(CFG, TEST_LOG_GROUP, stream_name),
        )

        current_time = Dates.now(tz"UTC")
        events = [
            LogEvent("First hey", current_time),
            LogEvent("Second hey", current_time + Second(1)),
            LogEvent("Third hey", current_time - Second(1)),
        ]

        setlevel!(LOGGER, "debug") do
            @test_log LOGGER "debug" "sorted" begin
                @test submit_logs(stream, events) == 3
            end
        end

        sleep(1)  # wait until AWS has injested the logs; this may or may not be enough
        response = logs(
            CFG,
            "GetLogEvents";
            logGroupName=TEST_LOG_GROUP,
            logStreamName=stream_name,
            startFromHead=true,
        )

        @test length(response["events"]) == 3
        messages = [event["message"] for event in response["events"]]
        @test messages == ["Third hey", "First hey", "Second hey"]
        delete_stream(CFG, TEST_LOG_GROUP, stream_name)
    end

    @testset "Rejected Logs" begin
        stream_name = new_stream("out_of_bounds")
        stream = CloudWatchLogStream(
            CFG,
            TEST_LOG_GROUP,
            create_stream(CFG, TEST_LOG_GROUP, stream_name),
        )

        current_time = Dates.now(tz"UTC")

        @test_warn LOGGER "retention policy" begin
            @test submit_log(stream, LogEvent("Way old", current_time - Day(8))) == 0
        end

        @test_warn LOGGER "days old" begin
            @test submit_log(stream, LogEvent("Too old", current_time - Day(15))) == 0
        end

        @test_warn LOGGER "hours in the future" begin
            @test submit_log(stream, LogEvent("Too new", current_time + Hour(3))) == 0
        end

        delete_stream(CFG, TEST_LOG_GROUP, stream_name)
    end

    @testset "Stream not found" begin
        @test_throws LOGGER StreamNotFoundException CloudWatchLogStream(
            CFG,
            TEST_LOG_GROUP,
            "made-up-stream-that-doesnt-exist",
        )
    end
end

@testset "CloudWatchLogHandler" begin
    @testset "Normal logging" begin
        start_time = CloudWatchLogs.unix_timestamp_ms()
        stream_name = new_stream("handler_type")

        handler = setlevel!(LOGGER, "debug") do
            @test_log LOGGER "debug" "initiated" begin
                handler = CloudWatchLogHandler(
                    CFG,
                    TEST_LOG_GROUP,
                    create_stream(CFG, TEST_LOG_GROUP, stream_name),
                    DefaultFormatter("{level} | {msg}"),
                )
                sleep(1)
            end

            handler
        end
        logger = Logger("CWLHLive.Normal"; propagate=false)
        push!(logger, handler)

        info(logger, "First log")
        warn(logger, "Second log")

        sleep(1)   # wait for the handler to submit the logs
        @test !isready(handler.channel)

        sleep(1)  # wait until AWS has injested the logs; this may or may not be enough
        response = logs(
            CFG,
            "GetLogEvents";
            logGroupName=TEST_LOG_GROUP,
            logStreamName=stream_name,
            startFromHead=true,
        )

        time_range = (start_time - 10):(CloudWatchLogs.unix_timestamp_ms() + 10)

        @test length(response["events"]) == 2
        messages = map(response["events"]) do event
            @test round(Int64, event["timestamp"]) in time_range
            event["message"]
        end

        @test messages == ["info | First log", "warn | Second log"]

        setlevel!(LOGGER, "debug") do
            # should cause the task to terminate with a debug message
            @test_log LOGGER "debug" "terminated normally" begin
                close(handler.channel)
                sleep(1)
            end
        end

        delete_stream(CFG, TEST_LOG_GROUP, stream_name)
    end

    @testset "Big logs" begin
        stream_name = new_stream("handler_big_logs")
        handler = CloudWatchLogHandler(
            CFG,
            TEST_LOG_GROUP,
            create_stream(CFG, TEST_LOG_GROUP, stream_name),
            DefaultFormatter("{msg}"),
        )
        logger = Logger("CWLHLive.Big"; propagate=false)
        push!(logger, handler)

        log_size = CloudWatchLogs.MAX_EVENT_SIZE - 26

        num_events = 26
        for (c, i) in zip('A':'Z', 1:num_events)
            info(logger, "$c" ^ log_size)
        end

        # wait for the logs to be submitted and for AWS to injest them
        sleep(10)
        response = logs(
            CFG,
            "GetLogEvents";
            logGroupName=TEST_LOG_GROUP,
            logStreamName=stream_name,
        )
        prev_token = ""
        num_events_injested = 0
        while prev_token != response["nextBackwardToken"]
            prev_token = response["nextBackwardToken"]
            @test length(response["events"]) <= 4
            num_events_injested += length(response["events"])

            response = logs(
                CFG,
                "GetLogEvents";
                logGroupName=TEST_LOG_GROUP,
                logStreamName=stream_name,
                nextToken=prev_token,
            )
        end
        @test num_events_injested == num_events
        delete_stream(CFG, TEST_LOG_GROUP, stream_name)
    end

    @testset "So many logs" begin
        start_time = CloudWatchLogs.unix_timestamp_ms()
        stream_name = new_stream("handler_so_many_logs")
        handler = CloudWatchLogHandler(
            CFG,
            TEST_LOG_GROUP,
            create_stream(CFG, TEST_LOG_GROUP, stream_name),
            DefaultFormatter("{msg}"),
        )
        logger = Logger("CWLHLive.SoMany"; propagate=false)
        push!(logger, handler)

        # not sure if this will actually go over the batch limit but we'll try
        max_num = CloudWatchLogs.MAX_BATCH_LENGTH * 2
        for i = 1:max_num
            info(logger, "$i")
        end

        # wait for the logs to be submitted
        for delay in CloudWatchLogs.PUTLOGEVENTS_DELAYS
            isready(handler.channel) || break
            sleep(delay)
        end

        sleep(1)  # wait until AWS has injested the logs; this may or may not be enough
        response = logs(
            CFG,
            "GetLogEvents";
            logGroupName=TEST_LOG_GROUP,
            logStreamName=stream_name,
            startFromHead=true,
            limit=1,
        )

        @test length(response["events"]) == 1
        event = response["events"][1]
        @test event["message"] == "1"

        sleep(5)  # wait until AWS has injested the logs; this may or may not be enough
        response = logs(
            CFG,
            "GetLogEvents";
            logGroupName=TEST_LOG_GROUP,
            logStreamName=stream_name,
            startFromHead=false,
            limit=1,
        )

        @test length(response["events"]) == 1
        event = response["events"][1]
        @test event["message"] == "$max_num"

        delete_stream(CFG, TEST_LOG_GROUP, stream_name)
    end
end

end


text:
function (::Type{QuadratureRule{2,RefTetrahedron}})(quad_type::Strang, order::Int)
    if order == 0 || order == 1
        points = [Vec{2}([1.0/3.0, 1.0/3.0])]
        weigths = [0.5]
    elseif order == 2
        points = [Vec{2}([1.0/6.0, 1.0/6.0]),
                  Vec{2}([1.0/6.0, 2.0/3.0]),
                  Vec{2}([2.0/3.0, 1.0/6.0])]
        weigths = [1/6,1/6,1/6]
    elseif order == 3
        points = [Vec{2}([0.659027622374092, 0.231933368553031]),
                  Vec{2}([0.659027622374092, 0.109039009072877]),
                  Vec{2}([0.231933368553031, 0.659027622374092]),
                  Vec{2}([0.231933368553031, 0.109039009072877]),
                  Vec{2}([0.109039009072877, 0.659027622374092]),
                  Vec{2}([0.109039009072877, 0.231933368553031])]
        weigths = 1/12*ones(Float64,6)
    elseif order == 4
        points = [Vec{2}([0.816847572980459, 0.091576213509771]),
                  Vec{2}([0.091576213509771, 0.816847572980459]),
                  Vec{2}([0.091576213509771, 0.091576213509771]),
                  Vec{2}([0.108103018168070, 0.445948490915965]),
                  Vec{2}([0.445948490915965, 0.108103018168070]),
                  Vec{2}([0.445948490915965, 0.445948490915965])]
        weigths = 1/2*ones(Float64,6)
        weigths[1:3] = 0.109951743655322*weigths[1:3]
        weigths[4:6] = 0.223381589678011*weigths[4:6]
    elseif order == 5
        points = [Vec{2}([0.33333333333333333, 0.33333333333333333]),
                  Vec{2}([0.79742698535308720, 0.10128650732345633]),
                  Vec{2}([0.10128650732345633, 0.79742698535308720]),
                  Vec{2}([0.10128650732345633, 0.10128650732345633]),
                  Vec{2}([0.05971587178976981, 0.47014206410511505]),
                  Vec{2}([0.47014206410511505, 0.05971587178976981]),
                  Vec{2}([0.47014206410511505, 0.47014206410511505])]
        weigths = 1/2*ones(Float64,7)
        weigths[1] = 0.22500000000000000*weigths[1]
        weigths[2:4] = 0.12593918054482717*weigths[2:4]
        weigths[5:7] = 0.13239415278850616*weigths[5:7]
    elseif order == 6
        points =[Vec{2}([0.873821971016996, 0.063089014491502]),
                 Vec{2}([0.063089014491502, 0.873821971016996]),
                 Vec{2}([0.063089014491502, 0.063089014491502]),
                 Vec{2}([0.501426509658179, 0.249286745170910]),
                 Vec{2}([0.249286745170910, 0.501426509658179]),
                 Vec{2}([0.249286745170910, 0.249286745170910]),
                 Vec{2}([0.636502499121399, 0.310352451033785]),
                 Vec{2}([0.636502499121399, 0.053145049844816]),
                 Vec{2}([0.310352451033785, 0.636502499121399]),
                 Vec{2}([0.310352451033785, 0.053145049844816]),
                 Vec{2}([0.053145049844816, 0.636502499121399]),
                 Vec{2}([0.053145049844816, 0.310352451033785])]
        weigths = 1/2*ones(Float64,12)
        weigths[1:3] = 0.050844906370207*weigths[1:3]
        weigths[4:6] = 0.116786275726379*weigths[4:6]
        weigths[7:12] = 0.082851075618374*weigths[7:12]
    else
        throw(ArgumentError("Strang rule of order $order not available"))
    end
    return QuadratureRule{2,RefTetrahedron,Float64}(weigths, points)
end