text:
using DiffEqBase, ParameterizedFunctions, DiffEqBiological, Markdown
#SDE Example Problems
export prob_sde_wave, prob_sde_linear, prob_sde_cubic, prob_sde_2Dlinear,
      prob_sde_lorenz, prob_sde_2Dlinear, prob_sde_additive,
      prob_sde_additivesystem, oval2ModelExample, prob_sde_bistable,
      prob_sde_bruss, prob_sde_oscilreact
#SDE Stratonovich Example Problems
export prob_sde_linear_stratonovich, prob_sde_2Dlinear_stratonovich

### SDE Examples

f_linear(u,p,t) = 1.01u
σ_linear(u,p,t) = 0.87u
linear_analytic(u0,p,t,W) = @.(u0*exp(0.63155t+0.87W))

@doc doc"""
```math
du_t = βudt + αudW_t
```
where β=1.01, α=0.87, and initial condtion u0=1/2, with solution

```math
u(u0,p,t,W_t)=u0\\exp((α-\\frac{β^2}{2})t+βW_t)
```

"""
prob_sde_linear = SDEProblem(SDEFunction(f_linear,σ_linear,
                             analytic=linear_analytic),σ_linear,1/2,(0.0,1.0))

linear_analytic_strat(u0,p,t,W) = @.(u0*exp(1.01t+0.87W))
prob_sde_linear_stratonovich = SDEProblem(SDEFunction(f_linear,σ_linear,
                               analytic=linear_analytic_strat),
                               σ_linear,1/2,(0.0,1.0))
f_linear_iip(du,u,p,t) = @.(du = 1.01*u)
σ_linear_iip(du,u,p,t) = @.(du = 0.87*u)
@doc doc"""
8 linear SDEs (as a 4x2 matrix):

```math
du_t = βudt + αudW_t
```
where β=1.01, α=0.87, and initial condtion u0=1/2 with solution

```math
u(u0,p,t,W_t)=u0\\exp((α-\\frac{β^2}{2})t+βW_t)
```
"""
prob_sde_2Dlinear = SDEProblem(SDEFunction(f_linear_iip,σ_linear_iip,
                               analytic=linear_analytic),
                               σ_linear_iip,ones(4,2)/2,(0.0,1.0))
prob_sde_2Dlinear_stratonovich = SDEProblem(SDEFunction(f_linear_iip,σ_linear_iip,
                               analytic=linear_analytic_strat),
                               σ_linear_iip,ones(4,2)/2,(0.0,1.0))

f_cubic(u,p,t) = -.25*u*(1-u^2)
σ_cubic(u,p,t) = .5*(1-u^2)
cubic_analytic(u0,p,t,W) = @. ((1+u0)*exp(W)+u0-1)/((1+u0)*exp(W)+1-u0)
ff_cubic = SDEFunction(f_cubic,σ_cubic,analytic = cubic_analytic)
@doc doc"""
```math
du_t = \\frac{1}{4}u(1-u^2)dt + \\frac{1}{2}(1-u^2)dW_t
```

and initial condtion u0=1/2, with solution

```math
u(u0,p,t,W_t)=\\frac{(1+u0)\\exp(W_t)+u0-1}{(1+u0)\\exp(W_t)+1-u0}
```
"""
prob_sde_cubic = SDEProblem(ff_cubic,σ_cubic,1/2,(0.0,1.0))

f_wave(u,p,t) = @. -0.01*sin(u)*cos(u)^3
σ_wave(u,p,t) = @. 0.1*cos(u)^2
wave_analytic(u0,p,t,W) = @. atan(0.1*W + tan(u0))
ff_wave = SDEFunction(f_wave,σ_wave,analytic=wave_analytic)
@doc doc"""
```math
du_t = -\\frac{1}{100}\sin(u)\cos^3(u)dt + \\frac{1}{10}\cos^{2}(u_t) dW_t
```

and initial condition `u0=1.0` with solution

```math
u(u0,p,t,W_t)=\\arctan(\\frac{W_t}{10} + \\tan(u0))
```
"""
prob_sde_wave = SDEProblem(ff_wave,σ_wave,1.0,(0.0,1.0))

f_additive(u,p,t) = @. p[2]/sqrt(1+t) - u/(2*(1+t))
σ_additive(u,p,t) = @. p[1]*p[2]/sqrt(1+t)
p = (0.1,0.05)
additive_analytic(u0,p,t,W) = @. u0/sqrt(1+t) + p[2]*(t+p[1]*W)/sqrt(1+t)
ff_additive = SDEFunction(f_additive,σ_additive,analytic=additive_analytic)
@doc doc"""
Additive noise problem

```math
u_t = (\\frac{β}{\\sqrt{1+t}}-\\frac{1}{2(1+t)}u_t)dt + \\frac{αβ}{\\sqrt{1+t}}dW_t
```

and initial condition u0=1.0 with α=0.1 and β=0.05, with solution

```math
u(u0,p,t,W_t)=\\frac{u0}{\\sqrt{1+t}} + \\frac{β(t+αW_t)}{\\sqrt{1+t}}
```
"""
prob_sde_additive = SDEProblem(ff_additive,σ_additive,1.0,(0.0,1.0),p)

f_additive_iip(du,u,p,t) = @.(du = p[2]/sqrt(1+t) - u/(2*(1+t)))
σ_additive_iip(du,u,p,t) = @.(du = p[1]*p[2]/sqrt(1+t))
ff_additive_iip = SDEFunction(f_additive_iip,σ_additive_iip,analytic=additive_analytic)
p = ([0.1;0.1;0.1;0.1],[0.5;0.25;0.125;0.1115])
@doc doc"""
A multiple dimension extension of `additiveSDEExample`

"""
prob_sde_additivesystem = SDEProblem(ff_additive_iip,σ_additive_iip,[1.;1.;1.;1.],(0.0,1.0),p)

f_lorenz = @ode_def_bare LorenzSDE begin
  dx = σ*(y-x)
  dy = x*(ρ-z) - y
  dz = x*y - β*z
end σ ρ β
σ_lorenz(du,u,p,t) = @.(du = 3.0)
@doc doc"""
Lorenz Attractor with additive noise

```math
\\begin{align}
dx &= σ*(y-x)dt + αdW_t \\\\
dy &= (x*(ρ-z) - y)dt + αdW_t \\\\
dz &= (x*y - β*z)dt + αdW_t \\\\
\\end{align}
```

with ``σ=10``, ``ρ=28``, ``β=8/3``, ``α=3.0`` and inital condition ``u0=[1;1;1]``.
"""
prob_sde_lorenz = SDEProblem(f_lorenz,σ_lorenz,ones(3),(0.0,10.0),(10.0,28.0,2.66))

f_nltest(u,p,t) = (1/3)*u^(1/3) + 6*u^(2/3)
σ_nltest(u,p,t) = u^(2/3)
analytic_nltest(u0,p,t,W) = (2t + 1 + W/3)^3
ff_nltest = SDEFunction(f_nltest,σ_nltest,analytic=analytic_nltest)
@doc doc"""
Runge–Kutta methods for numerical solution of stochastic differential equations
Tocino and Ardanuy
"""
prob_sde_nltest = SDEProblem(ff_nltest,σ_nltest,1.0,(0.0,10.0))

function oval2ModelExample(;largeFluctuations=false,useBigs=false,noiseLevel=1)
  #Parameters
  J1_200=3.
  J1_34=0.15
  J2_200=0.2
  J2_34=0.35
  J_2z=0.9
  J_O=0.918
  J_SO=0.5
  # J_ZEB=0.06
  J_ecad1=0.1
  J_ecad2=0.3
  J_ncad1=0.4
  J_ncad2=0.4
  J_ncad3 = 2
  J_snail0=0.6
  J_snail1=1.8
  J_zeb=3.0
  K1=1.0
  K2=1.0
  K3=1.0
  K4=1.0
  K5=1.0
  KTGF=20.
  Ks=100.
  # TGF0=0
  TGF_flg=0.
  Timescale=1000.
  dk_ZR1=0.5
  dk_ZR2=0.5
  dk_ZR3=0.5
  dk_ZR4=0.5
  dk_ZR5=0.5
  k0O=0.35
  k0_200=0.0002
  k0_34=0.001
  k0_snail=0.0005
  k0_zeb=0.003
  kO=1.2
  kOp=10.
  k_200=0.02
  k_34=0.019
  k_OT=1.1
  k_SNAIL=16.
  k_TGF=1.5
  k_ZEB=16.
  k_ecad0=5.
  k_ecad1=15.
  k_ecad2=5.
  k_ncad0=5.
  k_ncad1=2.
  k_ncad2=5.
  k_snail=0.05
  k_tgf=0.05
  k_zeb=0.06
  kdO=1.
  kd_200=0.035
  kd_34=0.035
  kd_SNAIL=1.6
  kd_SR1=0.9
  kd_TGF=0.9
  kd_ZEB=1.66
  kd_ecad=0.05
  kd_ncad=0.05
  kd_snail=0.09
  kd_tgf=0.1
  kd_tgfR=1.0
  kd_zeb=0.1
  kd_Op = 10.
  lamda1=0.5
  lamda2=0.5
  lamda3=0.5
  lamda4=0.5
  lamda5=0.5
  lamdas=0.5
  lamdatgfR=0.8
  nO=6.
  nSO=2.
  nzo=2.
  GE = 1.
  f = function (dy,y,p,t)
    # y(1) = snailt
    # y(2) = SNAIL
    # y(3) = miR34t
    # y(4) = SR1 # abundance of SNAIL/miR34 complex
    # y(5) = zebt
    # y(6) = ZEB
    # y(7) = miR200t
    # y(8) = ZR1 # abundance of ZEB/miR200 complex with i copies of miR200 bound on the sequence of ZEB1
    # y(9) = ZR2
    # y(10) = ZR3
    # y(11) = ZR4
    # y(12) = ZR5
    # y(13) = tgft
    # y(14) = TGF
    # y(15) = tgfR # abundance of TGF/miR200 complex
    # y(16) = Ecad
    # y(17) = Ncad
    # y(18) = Ovol2
    TGF0=.5(t>100)
    #ODEs
    dy[1]=k0_snail+k_snail*(((y[14]+TGF0)/J_snail0))^2/(1+(((y[14]+TGF0)/J_snail0))^2+(y[19]/J_SO)^nSO)/(1+y[2]/J_snail1)-kd_snail*(y[1]-y[4])-kd_SR1*y[4]
    dy[2]=k_SNAIL*(y[1]-y[4])-kd_SNAIL*y[2]
    dy[3]=k0_34+k_34/(1+((y[2]/J1_34))^2+((y[6]/J2_34))^2)-kd_34*(y[3]-y[4])-kd_SR1*y[4]+lamdas*kd_SR1*y[4]
    dy[4]=Timescale*(Ks*(y[1]-y[4])*(y[3]-y[4])-y[4])
    dy[5]=k0_zeb+k_zeb*((y[2]/J_zeb))^2/(1+((y[2]/J_zeb))^2+((y[19]/J_2z))^nO)-kd_zeb*(y[5]-(5*y[8]+10*y[9]+10*y[10]+5*y[11]+y[12]))-dk_ZR1*5*y[8]-dk_ZR2*10*y[9]-dk_ZR3*10*y[10]-dk_ZR4*5*y[11]-dk_ZR5*y[12]
    dy[6]=k_ZEB*(y[5]-(5*y[8]+10*y[9]+10*y[10]+5*y[11]+y[12]))-kd_ZEB*y[6]
    dy[7]=k0_200+k_200/(1+((y[2]/J1_200))^3+((y[6]/J2_200))^2)-kd_200*(y[7]-(5*y[8]+2*10*y[9]+3*10*y[10]+4*5*y[11]+5*y[12])-y[15])-dk_ZR1*5*y[8]-dk_ZR2*2*10*y[9]-dk_ZR3*3*10*y[10]-dk_ZR4*4*5*y[11]-dk_ZR5*5*y[12]+lamda1*dk_ZR1*5*y[8]+lamda2*dk_ZR2*2*10*y[9]+lamda3*dk_ZR3*3*10*y[10]+lamda4*dk_ZR4*4*5*y[11]+lamda5*dk_ZR5*5*y[12]-kd_tgfR*y[15]+lamdatgfR*kd_tgfR*y[15]
    dy[8]=Timescale*(K1*(y[7]-(5*y[8]+2*10*y[9]+3*10*y[10]+4*5*y[11]+5*y[12])-y[15])*(y[5]-(5*y[8]+10*y[9]+10*y[10]+5*y[11]+y[12]))-y[8])
    dy[9]=Timescale*(K2*(y[7]-(5*y[8]+2*10*y[9]+3*10*y[10]+4*5*y[11]+5*y[12])-y[15])*y[8]-y[9])
    dy[10]=Timescale*(K3*(y[7]-(5*y[8]+2*10*y[9]+3*10*y[10]+4*5*y[11]+5*y[12])-y[15])*y[9]-y[10])
    dy[11]=Timescale*(K4*(y[7]-(5*y[8]+2*10*y[9]+3*10*y[10]+4*5*y[11]+5*y[12])-y[15])*y[10]-y[11])
    dy[12]=Timescale*(K5*(y[7]-(5*y[8]+2*10*y[9]+3*10*y[10]+4*5*y[11]+5*y[12])-y[15])*y[11]-y[12])
    dy[13]=k_tgf-kd_tgf*(y[13]-y[15])-kd_tgfR*y[15]
    dy[14]=k_OT+k_TGF*(y[13]-y[15])-kd_TGF*y[14]
    dy[15]=Timescale*(TGF_flg+KTGF*(y[7]-(5*y[8]+2*10*y[9]+3*10*y[10]+4*5*y[11]+5*y[12])-y[15])*(y[13]-y[15])-y[15])
    dy[16]=GE*(k_ecad0+k_ecad1/(((y[2]/J_ecad1))^2+1)+k_ecad2/(((y[6]/J_ecad2))^2+1)-kd_ecad*y[16])
    dy[17]=k_ncad0+k_ncad1*(((y[2]/J_ncad1))^2)/(((y[2]/J_ncad1))^2+1)+k_ncad2*(((y[6]/J_ncad2))^2)/(((y[6]/J_ncad2)^2+1)*(1+y[19]/J_ncad3))-kd_ncad*y[17]
    dy[18]=k0O+kO/(1+((y[6]/J_O))^nzo)-kdO*y[18]
    dy[19]=kOp*y[18]-kd_Op*y[19]
  end

  σ1 = function (dσ,y,p,t)
    dσ[1] = noiseLevel*1.5y[1]
    dσ[18]= noiseLevel*6y[18]
  end

  σ2 = function (dσ,y,p,t)
    dσ[1] = 0.02y[1]
    dσ[16]= 0.02y[16]
    dσ[18]= 0.2y[18]
    dσ[17]= 0.02y[17]
  end

  if largeFluctuations
    σ = σ1
  else
    σ = σ2
  end

  u0 = [0.128483;1.256853;0.0030203;0.0027977;0.0101511;0.0422942;0.2391346;0.0008014;0.0001464;2.67e-05;4.8e-6;9e-7;0.0619917;1.2444292;0.0486676;199.9383546;137.4267984;1.5180203;1.5180203] #Fig 9B
  if useBigs
    u0 = big.(u0)
  end
  #u0 =  [0.1701;1.6758;0.0027;0.0025;0.0141;0.0811;0.1642;0.0009;0.0001;0.0000;0.0000;0.0000;0.0697;1.2586;0.0478;194.2496;140.0758;1.5407;1.5407] #Fig 9A
  SDEProblem(f,σ,u0,(0.0,500.0))
end

stiff_quad_f_ito(u,p,t) = -(p[1]+(p[2]^2)*u)*(1-u^2)
stiff_quad_f_strat(u,p,t) = -p[1]*(1-u^2)
stiff_quad_g(u,p,t) = p[2]*(1-u^2)

function stiff_quad_f_ito_analytic(u0,p,t,W)
    α = p[1]
    β = p[2]
    exp_tmp = exp(-2*α*t+2*β*W)
    tmp = 1 + u0
    (tmp*exp_tmp + u0 - 1)/(tmp*exp_tmp - u0 + 1)
end

ff_stiff_quad_ito = SDEFunction(stiff_quad_f_ito,stiff_quad_g,
                                analytic=stiff_quad_f_ito_analytic)

function stiff_quad_f_strat_analytic(u0,p,t,W)
    α = p[1]
    β = p[2]
    exp_tmp = exp(-2*α*t+2*β*W)
    tmp = 1 + u0
    (tmp*exp_tmp + u0 - 1)/(tmp*exp_tmp - u0 + 1)
end

ff_stiff_quad_strat = SDEFunction(stiff_quad_f_strat,stiff_quad_g,
                                analytic=stiff_quad_f_strat_analytic)

@doc doc"""
The composite Euler method for stiff stochastic
differential equations

Kevin Burrage, Tianhai Tian

And

S-ROCK: CHEBYSHEV METHODS FOR STIFF STOCHASTIC
DIFFERENTIAL EQUATIONS

ASSYR ABDULLE AND STEPHANE CIRILLI

Stiffness of Euler is determined by α+β²<1
Higher α or β is stiff, with α being deterministic stiffness and
β being noise stiffness (and grows by square).
"""
prob_sde_stiffquadito = SDEProblem(ff_stiff_quad_ito,stiff_quad_g,0.5,(0.0,3.0),(1.0,1.0))

@doc doc"""
The composite Euler method for stiff stochastic
differential equations

Kevin Burrage, Tianhai Tian

And

S-ROCK: CHEBYSHEV METHODS FOR STIFF STOCHASTIC
DIFFERENTIAL EQUATIONS

ASSYR ABDULLE AND STEPHANE CIRILLI

Stiffness of Euler is determined by α+β²<1
Higher α or β is stiff, with α being deterministic stiffness and
β being noise stiffness (and grows by square).
"""
prob_sde_stiffquadstrat = SDEProblem(ff_stiff_quad_strat,stiff_quad_g,0.5,(0.0,3.0),(1.0,1.0))

@doc doc"""
Stochastic Heat Equation with scalar multiplicative noise

S-ROCK: CHEBYSHEV METHODS FOR STIFF STOCHASTIC
DIFFERENTIAL EQUATIONS

ASSYR ABDULLE AND STEPHANE CIRILLI

Raising D or k increases stiffness
"""
function generate_stiff_stoch_heat(D=1,k=1;N = 100)
    A = full(Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1]))
    dx = 1/N
    A = D/(dx^2) * A
    function f(du,u,p,t)
        mul!(du,A,u)
    end
    #=
    function f(::Type{Val{:analytic}},u0,p,t,W)
        expm(A*t+W*I)*u0
    end
    =#
    function g(du,u,p,t)
        @. du = k*u
    end
    SDEProblem(f,g,ones(N),(0.0,3.0),noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM3)))
end

bistable_f(du,u,p,t) = du[1] = p[1] + p[2]*u[1]^4/(u[1]^4+11.9^4) - p[3]*u[1]
function bistable_g(du,u,p,t)
  du[1,1] = .1*sqrt(p[1] + p[2]*u[1]^4/(u[1]^4+11.9^4))
  du[1,2] = -.1*sqrt(p[3]*u[1])
end
p = (5.0,18.0,1.0)
"""
Bistable chemical reaction network with a semi-stable lower state.
"""
prob_sde_bistable = SDEProblem(bistable_f,bistable_g,[3.],(0.,300.),p,
                               noise_rate_prototype=zeros(1,2))

function bruss_f(du,u,p,t)
    du[1] = p[1] + p[2]*u[1]*u[1]*u[2] - p[3]*u[1]-p[4]*u[1]
    du[2] = p[3]*u[1] - p[2]*u[1]*u[1]*u[2]
end
function bruss_g(du,u,p,t)
  du[1,1] = 0.15*sqrt(p[1])
  du[1,2] = 0.15*sqrt(p[2]*u[1]*u[1]*u[2])
  du[1,3] = -0.15*sqrt(p[3]*u[1])
  du[1,4] = -0.15*sqrt(p[4]*u[1])
  du[2,1] = 0
  du[2,2] = -0.15*sqrt(p[2]*u[1]*u[1]*u[2])
  du[2,3] = 0.15*sqrt(p[3]*2.5*u[1])
  du[2,4] = 0
end
p = (1.0,1.0,2.5,1.0)
"""
Stochastic Brusselator
"""
prob_sde_bruss = SDEProblem(bruss_f,bruss_g,[3.,2.],(0.,100.),p,
                            noise_rate_prototype=zeros(2,4))

network = @reaction_network rnType  begin
    p1, (X,Y,Z) --> 0
    hill(X,p2,100.,-4), 0 --> Y
    hill(Y,p3,100.,-4), 0 --> Z
    hill(Z,p4,100.,-4), 0 --> X
    hill(X,p5,100.,6), 0 --> R
    hill(Y,p6,100.,4)*0.002, R --> 0
    p7, 0 --> S
    R*p8, S --> SP
    p9, SP + SP --> SP2
    p10, SP2 --> 0
end p1 p2 p3 p4 p5 p6 p7 p8 p9 p10
p = (0.01,3.0,3.0,4.5,2.,15.,20.0,0.005,0.01,0.05)

"""
An oscillatory chemical reaction system
"""
prob_sde_oscilreact = SDEProblem(network, [200.,60.,120.,100.,50.,50.,50.], (0.,4000.),p)


text:
"""
    Broadcasting()

Broadcast inner reducing function over elements in the input.  Roughly
speaking, it transforms the inner reducing function `op` to `op′(a, b)
= op.(a, b)`.  However, it has a better memory usage and better
initial value handling.

If the input is an array, the array created at the first iteration is
reused if it can hold the element types of subsequent iterations.
Otherwise, the array type is widen as needed.

If `init` passed to the fold function is a lazy "initializer" object
such as [`OnInit`](@ref), it is initialized independently for each
item in the first input array.  This makes using `Broadcasting` for
(possibly) in-place functions safe.

!!! compat "Transducers.jl 0.4.32"

    New in version 0.4.32.

!!! note

    `Broadcasting` transducer is not supported in Julia 1.0.

# Examples
```jldoctest Broadcasting
julia> using Transducers

julia> foldl(+, Broadcasting(), [[1, 2], [3, 4], 5])
2-element Vector{Int64}:
  9
 11

julia> foldl(+, Broadcasting(), [(0,), [1], [2.0], [3 + 0im]])
1-element Vector{ComplexF64}:
 6.0 + 0.0im

julia> foldl(
           *,
           [[[1], [10, 100]], [[2], [20, 200]], [[3], [30, 300]]] |>
               Broadcasting() |> Broadcasting(),
       )
2-element Vector{Vector{Int64}}:
 [6]
 [6000, 6000000]
```

When processing nested data structure (e.g., array-of-arrays) and
mutating accumulator in-place, be careful with sharing the
accumulators with each processing of items in the input.  For example,
this is a bad example:

```jldoctest Broadcasting
julia> add!(a, b) = a .+= b;

julia> foldl(add!, Broadcasting(), [[[1], [2, 3]], [[4, 5], 6]];
             init = Ref([0, 0]))
2-element Vector{Vector{Int64}}:
 [13, 15]
 [13, 15]

julia> ans[1] === ans[2]  # they are the same object
true
```

Use `OnInit` to initialize a new array with each item in the input:

```jldoctest Broadcasting
julia> foldl(add!, Broadcasting(), [[[1], [2, 3]], [[4, 5], 6]];
             init = OnInit(() -> [0, 0]))
2-element Vector{Vector{Int64}}:
 [5, 6]
 [8, 9]

julia> ans == [
           add!(add!([0, 0], [1]), [4, 5]),
           add!(add!([0, 0], [2, 3]), 6),
       ]
true
```
"""
struct Broadcasting <: Transducer end

OutputSize(::Type{<:Broadcasting}) = SizeStable()

struct LazyInit{I}
    init::I
end

# Defer calling inner `start`. This is for broadcasting `init` to the
# shape of the first item:
start(::R_{Broadcasting}, init) = LazyInit(init)

@inline next(rf::R_{Broadcasting}, acc::LazyInit, xs) =
    next.(inner(rf), start.(inner(rf), acc.init), xs)
@inline function next(rf::R_{Broadcasting}, accs, xs)
    @! accs .= next.(inner(rf), accs, xs)
    return accs
end
@inline (rf::R_{Broadcasting})(accs, xs) = next(rf, accs, xs)

complete(rf::R_{Broadcasting}, ::LazyInit{typeof(DefaultInit)}) =
    DefaultInitOf{typeof(rf)}()

# Is this the right thing to do?
@inline complete(rf::R_{Broadcasting}, acc::LazyInit) = start.(inner(rf), acc.init)

@inline function complete(rf::R_{Broadcasting}, accs)
    @! accs .= complete.(inner(rf), accs)
    if any(r -> r isa DefaultInitOf, accs)  # is `any` too strict?
        return DefaultInitOf{typeof(rf)}()
    end
    return accs
end

combine(::R_{Broadcasting}, ::LazyInit, accs) = accs
combine(::R_{Broadcasting}, accs, ::LazyInit) = accs
combine(::R_{Broadcasting}, accs::LazyInit, ::LazyInit) = accs
@inline function combine(rf::R_{Broadcasting}, lefts, rights)
    @! lefts .= combine.(inner(rf), lefts, rights)
    return lefts
end

ext: jl


text:
# Ground truth was generated using the bw.SJ() R function, with nb=10000 and tol=1e-3 unless further specified.
# Some differences are expected because of how the double sums are approximated.

@testset "fewpoints" begin
	@test isnan(bwsj([0.0]))
	@test isnan(bwsj([1.0]))
	@test isnan(bwsj(ones(100)))
	@test bwsj([1.0,2.0])≈0.1163591 rtol=0.1
end

@testset "range" begin
	# R: nb=10000000, tol=1e-6
	@test bwsj(range(0,stop=1,length= 96); rtol=1e-1)≈0.1109613 rtol=1e-1
	@test bwsj(range(0,stop=1,length= 97); rtol=1e-1)≈0.1105281 rtol=1e-1
	@test bwsj(range(0,stop=1,length= 98); rtol=1e-1)≈0.1101011 rtol=1e-1
	@test bwsj(range(0,stop=1,length= 99); rtol=1e-1)≈0.1096801 rtol=1e-1
	@test bwsj(range(0,stop=1,length=100); rtol=1e-1)≈0.1092651 rtol=1e-1
	@test bwsj(range(0,stop=1,length=101); rtol=1e-1)≈0.1088558 rtol=1e-1
	@test bwsj(range(0,stop=1,length=102); rtol=1e-1)≈0.1084521 rtol=1e-1
	@test bwsj(range(0,stop=1,length=103); rtol=1e-1)≈0.1080540 rtol=1e-1
	@test bwsj(range(0,stop=1,length=104); rtol=1e-1)≈0.1076612 rtol=1e-1
	@test bwsj(range(0,stop=1,length=105); rtol=1e-1)≈0.1072737 rtol=1e-1

	@test bwsj(range(0,stop=1,length= 96); rtol=1e-3)≈0.1109613 rtol=1e-3
	@test bwsj(range(0,stop=1,length= 97); rtol=1e-3)≈0.1105281 rtol=1e-3
	@test bwsj(range(0,stop=1,length= 98); rtol=1e-3)≈0.1101011 rtol=1e-3
	@test bwsj(range(0,stop=1,length= 99); rtol=1e-3)≈0.1096801 rtol=1e-3
	@test bwsj(range(0,stop=1,length=100); rtol=1e-3)≈0.1092651 rtol=1e-3
	@test bwsj(range(0,stop=1,length=101); rtol=1e-3)≈0.1088558 rtol=1e-3
	@test bwsj(range(0,stop=1,length=102); rtol=1e-3)≈0.1084521 rtol=1e-3
	@test bwsj(range(0,stop=1,length=103); rtol=1e-3)≈0.1080540 rtol=1e-3
	@test bwsj(range(0,stop=1,length=104); rtol=1e-3)≈0.1076612 rtol=1e-3
	@test bwsj(range(0,stop=1,length=105); rtol=1e-3)≈0.1072737 rtol=1e-3
end

@testset "reference" begin
	@test bwsj_reference(range(0,stop=1,length=100); rtol=1e-1)≈0.1092651 rtol=1e-1
end

@testset "scale" begin
	@test bwsj(range(0,stop=1e3,  length=101))≈0.1088518*1e3   rtol=0.1
	@test bwsj(range(0,stop=1e-3, length=101))≈0.1088518*1e-3  rtol=0.1
	@test bwsj(range(0,stop=1e24, length=101))≈0.1088518*1e24  rtol=0.1
	@test bwsj(range(0,stop=1e-24,length=101))≈0.1088518*1e-24 rtol=0.1
end

@testset "randn" begin
	rng = StableRNG(1234)
	# @info "A"
	X = randn(rng, 20)
	@test bwsj(X)≈0.6189722 rtol=0.1

	# @info "B"
	X = randn(rng, 100)
	@test bwsj(X)≈0.4314236 rtol=0.1

	# @info "C"
	X = randn(rng, 1000)
	@test bwsj(X)≈0.2823561 rtol=0.1
end

@testset "quantile" begin
	N = 100
	X = quantile.(Normal(), range(1/2N,stop=1-1/2N,length=N))

	@test bwsj(sign.(X).*abs.(X).^(2.0.^-4))≈0.05301544  rtol=0.1
	@test bwsj(sign.(X).*abs.(X).^(2.0.^-3))≈0.07877113  rtol=0.1
	@test bwsj(sign.(X).*abs.(X).^(2.0.^-2))≈0.1279672   rtol=0.1
	@test bwsj(sign.(X).*abs.(X).^(2.0.^-1))≈0.2092015   rtol=0.1
	@test bwsj(sign.(X).*abs.(X).^(2.0.^0 ))≈0.4747508   rtol=0.1
	@test bwsj(sign.(X).*abs.(X).^(2.0.^1 ))≈0.1394881   rtol=0.1
	@test bwsj(sign.(X).*abs.(X).^(2.0.^2 ))≈0.03132362  rtol=0.1
	@test bwsj(sign.(X).*abs.(X).^(2.0.^3 ))≈0.004705298 rtol=0.1 # NB: run with nb=10000000 to get accurate answer

	@test bwsj(sign.(X).*abs.(X).^(2.0.^6))≈1.1579829602522003e-12 rtol=0.1 # NB: ground truth generated with bwsj_reference(; rtol=1e-24)
end

@testset "outliers" begin
	# ground truth generated with bwsj_reference(; rtol=1e-24)
	rng = StableRNG(2017)
	X = rand(rng, 100)

	Y1 = 1e9 .+ randn(rng,1)
	@test bwsj(vcat(X,Y1))≈0.11706501880745836 rtol=0.1

	Y2 = 1e9 .+ randn(rng,10)
	@test bwsj(vcat(X,Y2))≈0.1197213803664104 rtol=0.1

	Y3 = 1e9.*sign.(randn(rng,10)) .+ randn(rng,10)
	@test bwsj(vcat(X,Y3))≈0.12247926574875453 rtol=0.1
end

@testset "repeats" begin
    X = vcat(zeros(10),ones(10),2*ones(10),4*ones(10))
    @test bwsj(X;leafsize=10)≈0.1246847 rtol=1
    @test bwsj(X;leafsize=5 )≈0.1246847 rtol=1
end

@testset "ϕbounds" begin
	for stepsize in (0.1, 0.5, 1.0, 2.5, 5.0)
		a = 0:stepsize:5
		b = a.+stepsize
		x = range(0,stop=stepsize,length=11)[2:end-1]'

		B = ϕ4bounds.(a,b,a.+x)
		lb,ub = first.(B),last.(B)
		@test all(lb.<=ϕ4.(a.+x).<=ub)

		B = ϕ6bounds.(a,b,a.+x)
		lb,ub = first.(B),last.(B)
		@test all(lb.<=ϕ6.(a.+x).<=ub)
	end
end


text:
function enso_classification(year)
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
module Stdlib
  #import Base: include
  export include

  function include(lib::String, file::String)
    if file == "*"
      # Load all files from lib
      if ispath(@__DIR__, "stdlib", lib, "includes.jl")
        # Use includes.jl if available
        Base.include(joinpath(@__DIR__, "stdlib", lib, "includes.jl"))
      else
        # If no includes.jl, load in any order
        for f in filter(x->endswith(x,".jl"), readdir(joinpath(@__DIR__, "stdlib", lib)))
          path = joinpath(@__DIR__, "stdlib", lib, f)
          Base.include(path)
        end
      end
    else
      # Load specified file from lib
      path = joinpath(@__DIR__, "stdlib", lib, file)
      @assert ispath(path) "Failed to locate $lib/$file in Stdlib"
      Base.include(path)
    end
  end
  function include_lib(lib::String)
    if lib == "*"
      # TODO: Load all libraries and files
    else
      # TODO: Load all files from specified lib
    end
  end
end


text:
let
    arr = ones(2,2)
    pad_matrix(arr,(2,2),0)
end
let
    a = [1.0,2,3]
    b = [-4.0,5,6]
    @test array_isapprox(cross(a,b),cross_product_operator(a)*b)
end
let
    n = 4
    for i in 1:n
        v = one_hot(n,i)
        @test v[i] == 1
        v = one_hot(i,n)
        @test v[i] == 1
        for j in 1:n
            if i != j
                @test v[j] == 0
            end
        end
    end
    @test_throws AssertionError one_hot(-1,3)
    @test_throws AssertionError one_hot(1,-3)
    @test_throws AssertionError one_hot(-1,-3)
end
let
    vec = [1,2,3,4]
    for i in 1:4
        @test wrap_idx(length(vec),i) == i
        @test wrap_get(vec,i)[1] == i
    end
    @test wrap_idx(length(vec),5) == 1
    @test wrap_get(vec,5)[1] == 1
    @test wrap_idx(length(vec),9) == 1
    @test wrap_idx(length(vec),0) == 4
    @test wrap_idx(length(vec),-4) == 4
    @test wrap_idx(length(vec),-5) == 3
end
let
    a = [
        1 2;
        3 4
        ]
    @test wrap_get(a,[1,1]) == 1
    @test wrap_get(a,[3,3]) == 1
    @test wrap_get(a,[4,4]) == 4
    @test wrap_get(a,[0,-1]) == 3
end


text:
using BlackBoxOptimizationBenchmarking, CMAEvolutionStrategy, Evolutionary, NaturalES
const BBOB = BlackBoxOptimizationBenchmarking
include(joinpath(pathof(BBOB), "..", "..", "scripts", "optimizers_interface.jl"))
pinit(D) = 10*rand(D).-5
c = opt -> Chain(opt, NelderMead(), 0.9)
function BBOB.optimize(m::Chain,f,D,run_length)
    rl1 = round(Int,m.p*run_length)
    rl2 = run_length - rl1
    mfit = BBOB.optimize(m.first,f,D,run_length)
    xinit = BBOB.minimizer(mfit)
    mfit = BBOB.optimize(m.second,f,D,run_length,xinit)
end

struct CMAES end
function BBOB.optimize(::CMAES, f, D, run_length)
    CMAEvolutionStrategy.minimize(f, pinit(D), 3., verbosity = 0, maxfevals = run_length)
end
BBOB.minimum(o::CMAEvolutionStrategy.Optimizer) = fbest(o)
BBOB.minimizer(o::CMAEvolutionStrategy.Optimizer) = xbest(o)
struct Evo end
function BBOB.optimize(o::Evo, f, D, run_length)
    lambda = CMAEvolutionStrategy.default_popsize(D)
    Evolutionary.optimize(f, pinit(D), Evolutionary.CMAES(λ = lambda), Evolutionary.Options(iterations = ceil(Int, run_length/lambda)))
end
BBOB.minimum(o::Evolutionary.EvolutionaryOptimizationResults) = o.minimum
BBOB.minimizer(o::Evolutionary.EvolutionaryOptimizationResults) = o.minimizer
struct NES{T} end
function BBOB.optimize(o::NES{T}, f, D, run_length) where T
    lambda = CMAEvolutionStrategy.default_popsize(D)
    NaturalES.optimize(f, pinit(D), 1., T, samples = lambda,
                       iterations = ceil(Int, run_length/lambda))
end
BBOB.minimum(o::NamedTuple{(:sol, :cost),Tuple{Array{Float64,1},Float64}}) = o.cost
BBOB.minimizer(o::NamedTuple{(:sol, :cost),Tuple{Array{Float64,1},Float64}}) = o.sol

ms = [
      c(BlackBoxOptimMethod(:adaptive_de_rand_1_bin_radiuslimited)),
      c(BlackBoxOptimMethod(:xnes)),
      c(NES{xNES}()),
      c(NES{sNES}()),
      c(Evo()),
      c(PyCMA()),
      c(CMAES())
      ]
D = 24
lengths = round.(Int,range(1_000, stop=80_000, length=5))
res = BBOB.benchmark(ms, 1:length(enumerate(BBOBFunction)), lengths, 2, D, 1e-6)

using Statistics, PGFPlotsX
mres = reshape(mean(res[1], dims = 2), 7, :)
mnames = ["BBO.default", "BBO.xnes", "NES.xNES", "NES.sNES", "Evo.CMAES", "PyCMA", "CMA-ES"]
plot1 = @pgf Axis({legend_entries = mnames, legend_pos = "outer north east", xticklabels = lengths, xtick = 1:length(lengths), font = raw"\small", xlabel = "run length", ylabel = "success rate", title = "All functions, D: 24"},
          [PlotInc(Coordinates(1:5, mres[i, :])) for i in 1:7]...)
pgfsave(joinpath(@__DIR__, "..", "bbob24.png"), plot1)
plot2 = @pgf Axis({ytick = 1:7, yticklabels = mnames, xlabel = "relative run time"},
          Plot({xbar}, Table(reshape(mean(res[4], dims = 2), :), 1:7)))
pgfsave(joinpath(@__DIR__, "..", "bbob24rt.png"), plot2)

max_stars_repo_name: jbrea/CMAEvolutionStrategy

max_issues_repo_name: jbrea/CMAEvolutionStrategy


text:
using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using StatsBase
using LinearAlgebra
using GLM
using Test
using Random
using BenchmarkTools
using DelimitedFiles
using BGEN
using VCFTools

# write your own tests here
include("utilities_test.jl")
include("L0_reg_test.jl")
include("multivariate_test.jl")
include("cv_iht_test.jl")
include("wrapper_test.jl")


text:
module NER
using Main, JLD, ArgParse, Knet

# opt: adam (beta1=0.9, beta2=0.999,epsilon=1e-08) batch=32, lr=0.001, gclip=5
# hidden: 128; no recurrent; biyofiz'de cpu'da 1 trn epoch 18 thread ile 258 sn. => 0.0371 char error
# hidden: 512; output layer'da recurrent connection; 1trn epoch: gpu'da 217 sn.; dev cerr (char error): 0.0256

# Knet speed: hidden:128, batch:32, non-recurr out, --fast 154 sec/epoch
# Knet speed: hidden:512, batch:32, non-recurr out, --fast 208 sec/epoch

isdefined(:xtrn) || (@load "ner.jld")

function main(args=ARGS)

    # Set training parameters:

    s = ArgParseSettings()
    @add_arg_table s begin
        ("--hidden"; default=128; arg_type=Int)
        ("--lr"; default=0.001; arg_type=Float64) # WARNING: adam not implemented yet
        ("--batchsize"; default=32; arg_type=Int) # WARNING: I skip the leftovers
        ("--epochs"; default=1000; arg_type=Int)
        ("--gclip"; default=5.0; arg_type=Float64; help="model-wide gradient clip")
        ("--winit"; default="Gaussian(0,0.1)"; help="weight initialization")
        ("--fast"; action=:store_true; help="skip norm and loss calculations.")
        ("--devfortrn"; action=:store_true; help="use dev for training for timing experiments.")
        ("--gcheck"; default=0; arg_type=Int; help="gradient checking")
        ("--lossreport"; default=0; arg_type=Int; help="report training loss every n tokens")
        ("--seed"; default=42; arg_type=Int; help="random seed")
    end
    isa(args, AbstractString) && (args=split(args))
    opts = parse_args(args, s)
    println(opts)
    for (k,v) in opts; @eval ($(symbol(k))=$v); end
    winit1 = eval(parse(winit))
    seed > 0 && setseed(seed)

    # Load data: (should we shuffle?)

    @show map(size, (xtrn, ytrn, xdev, ydev))
    global trn = TagData(xtrn, ytrn; batchsize=batchsize, dense=true)
    global dev = TagData(xdev, ydev; batchsize=batchsize, dense=true)
    @show maxtoken(xtrn)
    @show nclass = maxtoken(ytrn)
    flush(STDOUT)

    # Construct tagger model

    # initialization ek:
    # output layer'da sofmax kullanirsak:
    # W=lasagne.init.GlorotUniform(gain=1.0) # aka Xavier initialization.
    # b=lasagne.init.Constant(0.)

    # output layer'da recurrent softmax kullanirsak:
    # W_in_to_hid=lasagne.init. GlorotUniform(gain=1.0) # aka Xavier initialization.
    # W_hid_to_hid=Identity(),
    # b=lasagne.init.Constant(0.)
    # hid_init=lasagne.init.Constant(0.)

    global model = Tagger(peeplstm, peeplstm, pred; nclass=nclass, hidden=hidden, winit=winit1)
    setopt!(model; lr=lr)
    losscnt = (fast ? nothing : zeros(2))
    maxnorm = (fast ? nothing : zeros(2))
    history = Any[]
    epoch=trnprp=devprp=deverr=0

    for epoch=1:epochs
        @show epoch
        devfortrn ?
        train(model, dev, softloss; gclip=gclip, losscnt=losscnt, maxnorm=maxnorm, lossreport=lossreport) :
        train(model, trn, softloss; gclip=gclip, losscnt=losscnt, maxnorm=maxnorm, lossreport=lossreport)
        gcheck > 0 && gradcheck(model, trn, softloss; gcheck=gcheck)
        @show trnprp = exp(losscnt[1]/losscnt[2])
        @show devprp = exp(test(model, dev, softloss))
        @show deverr = test(model, dev, zeroone)
        push!(history, deverr)
        if length(history) > 5 && history[end] > history[end-5]
            @show lr /= 2
            setopt!(model; lr=lr)
        end
        flush(STDOUT)
    end

    return (epoch,trnprp,devprp,deverr)
end

# Construct peephole lstm model from http://arxiv.org/pdf/1308.0850.pdf pp.5
# TODO: faster lstm possible if I do it with a single wdot
# add2,wdot,bias defined in Knet/src/op/compound.jl

@knet function peeplstm(x; hidden=0, o...)
    input  = add3(x,h,c; o..., f=sigm, out=hidden)
    forget = add3(x,h,c; o..., f=sigm, out=hidden)
    newmem = add2(x,h;   o..., f=tanh, out=hidden)
    ig = mul(input,newmem)
    fc = mul(c,forget)
    c  = add(ig,fc)
    output = add3(x,h,c; o..., f=sigm, out=hidden)
    tc = tanh(c)
    h  = mul(output,tc)
end

@knet function add3(x1, x2, x3; f=sigm, o...)
    y1 = wdot(x1; o...)
    y2 = wdot(x2; o...)
    y3 = wdot(x3; o...)
    z1 = add(y2,y1)
    z2 = add(z1,y3)
    z3 = bias(z2; o...)
    ou = f(z3; o...)
end

@knet function pred(x,y; nclass=0, o...)
    z = add2(x,y; out=nclass, f=soft, winit=Xavier())
end

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

end # module

### DEAD CODE:

## We should use kfun's for constructor, not Nets
# fnet = Net(peeplstm; out=hidden, winit=winit)
# bnet = Net(peeplstm; out=hidden, winit=winit)
# pnet = Net(add2; ninputs=2, out=nclass, f=soft, winit=Xavier())
# model = Tagger(fnet, bnet, pnet)


text:
module GAIO

using LinearAlgebra
using SparseArrays
using StaticArrays
# using GLFW
# using ModernGL
using GeometryBasics
using WGLMakie
using LightGraphs
using ForwardDiff
using Arpack
using Base.Threads

using GLMakie

export Box

export AbstractBoxPartition, BoxPartition, TreePartition
export dimension, depth, key_to_box, point_to_key, tree_search

export BoxSet
export boxset_empty, subdivide, subdivide!

export BoxFun

export TransferOperator
export strongly_connected_components, matrix, eigs

export BoxMap, PointDiscretizedMap
export SampledBoxMap, AdaptiveBoxMap
export boxmap

export map_boxes, map_boxes_new

export rk4, rk4_flow_map

export relative_attractor, unstable_set!, chain_recurrent_set, cover_roots, finite_time_lyapunov_exponents

export plot

# ENV["JULIA_DEBUG"] = all

include("box.jl")

abstract type AbstractBoxPartition{B <: Box} end

include("partition_regular.jl")
include("partition_tree.jl")
include("boxset.jl")
include("boxmap.jl")
include("boxfun.jl")
include("transfer_operator.jl")
include("algorithms.jl")
include("plot.jl")

end # module


text:
using Mimi
using OptiMimi
using Test

@defcomp Simple begin
    region = Index()

    xx = Parameter(index=[time, region])
end

m = Model()
set_dimension!(m, :region, collect(1:2))
set_dimension!(m, :time, collect(1:3))

simple = add_comp!(m, Simple)
xx = reshape(repeat([1, 2], inner=3), (3, 2))
simple[:xx] = xx

gen(rr) = m.md.external_params[:xx].values[1, rr]
hall = hallsingle(m, :Simple, :xx, gen, [:time])

@test reshape(hall.f, (3, 2)) == m.md.external_params[:xx].values[:, :]


text:
#oneoff tasks
function oneOffCategories(data::NCCSData)
  metadf::DataFrame = deepcopy(data.df[data.df[:fisyr] .== 2013, [:category, :categoryname, :naics2name, :naics2]])
  dropmissing!(metadf)
  outdf::DataFrame = by(metadf, [:category, :categoryname, :naics2name, :naics2],
    [:category, :categoryname, :naics2name, :naics2] =>
      x -> (n=length(x.category)))

  CSV.write("$WORKING_PATH\\categoryMapping.csv",outdf)

  return nothing
end

function oneOffTasks(data::NCCSData)
  println(data.df[1:10, :datefiscal])

  return nothing
end

#conducts main data processing and analysis
#results are cached in the serial file. Useful to avoid redoing computation
function constructPrimaryOutput(; refreshFilter::Bool = true,
    primaryOutputName::String = PRIMARY_OUTPUT_NAME,
    writePrimaryCSV::Bool = false,
    dataPath::String = DATA_PATH,
    refreshoutsidedata::Bool = false,
    filternoncentral::Bool = true,
    refreshopen990data=true)::NCCSData

  local data::NCCSData
  local outPath = "$dataPath\\$(primaryOutputName)_long.csv"
  if refreshFilter
    data = deserializeNCCSData()
    data = processAll!(data,
      refreshoutsidedata=refreshoutsidedata,
      filternoncentral=filternoncentral,
      refreshopen990data=refreshopen990data)
    serializeNCCSData(data, dataName = primaryOutputName)
  else
    data = deserializeNCCSData(dataName = primaryOutputName)
  end

  #for debugging
  #uCSV.write(outPath, data.df[400_000:450_000,:])

  writePrimaryCSV && uCSV.write(outPath, data.df) #write to a csv if desired

  return data

end

function writedata(df::DataFrame; rows::AbstractRange = 1:size(df,1),
  cols::Vector{Symbol}=names(df),
  testname::String = "dataout", workingpath::String = WORKING_PATH)

  path::String = "$workingpath\\dataout.csv"
  CSV.write(path, df[rows,cols])
end

#launching point for the descriptive output
function reportandanalyze(;refreshFilter::Bool = true,
    refreshDescriptiveFigures::Bool = true,
    refreshDescriptiveTables::Bool = true, writePrimaryCSV::Bool = false,
    doOneOffTasks::Bool= false,
    refreshβ::Bool = true,
    refreshoutsidedata::Bool = false,
    refreshpersistence::Bool = false,
    refreshtopinstitutions::Bool = false,
    refreshopen990data::Bool = true,
    refreshsensitivity::Bool = true
    )

    @time data::NCCSData = constructPrimaryOutput(refreshFilter=refreshFilter,
      writePrimaryCSV=writePrimaryCSV,
      refreshoutsidedata=refreshoutsidedata,
      refreshopen990data=refreshopen990data)

    #NOTE: Begin debug code
    sort!(data.df, [:ein, :fisyr])

    #=writedata(data.df[data.df[:fisyr] .≥ 2009,[:ein, :fisyr, :lreturn,
      :Llcontributions, :lcontributions, :contributions,
      :Lpublicsupport, :publicsupport, :publicsupport_avg,
      :Llprogramexpense, :lprogramexpense, :programexpense,
      :programexpense_avg, :lprogramexpense_avg,
      :LpublicsupportXLlcont,  :LlprogramexpenseXLlcont, :publicsupport_avgXLlcont,
      :lprogramexpense_avgXLlcont
      ]])=#

    #end debug code

    refreshDescriptiveFigures && makeDescriptiveFigures(data)
    refreshDescriptiveTables && makeDescriptiveTables(data)
    refreshtopinstitutions && getTopInstitutions(data, sorted=true) #sorted=true ok after descriptive tables
    refreshβ && processβ(data)
    refreshpersistence && analyzepersistence(data)
    refreshsensitivity && olssensitivitytables(data)

    doOneOffTasks && oneOffTasks(data)

    return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Random
# Defining function constants
const BIG_NUMBER = 100
const CONVERG_ITER = 10
const CONVERG_CUTOFF = 0.005  # cosine distance
const TOTAL_INIT_CONDITIONS = 5

function evaluateStability(Wall, Hall, totalProcesses, totalReplicates)

    Wall = map(x->max(x, EPS), Wall)
    Hall = map(x->max(x, EPS), Hall)

    Random.seed!(123456)

    # Clustering mutational processes using custom clustering procedure
    minClusterDist = BIG_NUMBER
    totalIter = div(size(Wall, 2), totalProcesses)
    idx = zeros(size(Hall, 1))
    clusterCompactness = zeros(totalProcesses, totalIter)
    iStartigDataSet = 1 : totalProcesses : size(Wall,2)
    iStartingDataSet = iStartigDataSet[randperm(totalIter)]

    local centroids, centroidsFinal, clusterCompactnessFinal, idxFinal, centroidStd
    local exposure, exposureStd, idxS, processStab, processStabAvg

    for iInitData = 1 : min(TOTAL_INIT_CONDITIONS, totalIter) % size(Wall, 2)
        iStartingData = iStartingDataSet[iInitData]
        iEnd = iStartingData + totalProcesses - 1
        centroids = Wall[:, iStartingData:iEnd]

        # this is random process
        centroidsTest = rand(Float64, size(centroids))
        countIRep = 0

        for iRep = 1 : totalReplicates
            allDist = pairwise(CosineDist(), hcat(centroids, Wall), dims=2)
            centroidDist = allDist[ 1:size(centroids, 2), (size(centroids, 2)+1):size(allDist, 2) ]'

            jRange = randperm(totalProcesses)
            for jIndex = 1 : totalProcesses
                j = jRange[jIndex]
                for i = 1 : totalProcesses : size(Wall, 2)
                    iRange = i: (i + totalProcesses - 1)
                    Ind = argmin( centroidDist[iRange, j] )[1]
                    centroidDist[iRange[Ind], :] .= BIG_NUMBER
                    idx[iRange[Ind]] = j
                end
            end

            maxDistToNewCentroids = 0
            for i = 1 : size(centroids, 2)
                centroids[:, i] = mean( Wall[:, idx .== i], dims=2 )
                maxDistToNewCentroids = max(maxDistToNewCentroids, maximum(pairwise(CosineDist(), hcat(centroids[:, i], centroidsTest[:, i]), dims=2)))
            end

            if maxDistToNewCentroids < CONVERG_CUTOFF
                countIRep =  countIRep + 1
            else
                countIRep = 0
                centroidsTest = centroids
            end

            if countIRep == CONVERG_ITER
                break
            end
        end

        for i = 1 : size(centroids, 2)
            clusterDist = pairwise(CosineDist(), hcat(centroids[:,i], Wall[:, idx.==i]), dims=2)
            clusterCompactness[i, :] = clusterDist[1, 2:size(clusterDist, 2)]
        end

        meanCompactness = mean(clusterCompactness)
        if minClusterDist > meanCompactness
            centroidsFinal = centroids
            idxFinal = copy(idx)
            clusterCompactnessFinal = clusterCompactness
            minClusterDist = meanCompactness
        end
    end

    centroids = centroidsFinal'
    idx = copy(idxFinal)
    clusterCompactness = clusterCompactnessFinal

    centDist = mean( clusterCompactness, dims=2 )[:,1]
    centDistInd = sortperm(centDist)  # centDistSorted not used
    clusterCompactness = clusterCompactness[centDistInd,:]

    centroids = centroids[centDistInd, :]
    idxNew = zeros(size(idx));

    for i = 1 : totalProcesses
        idxNew[idx .== centDistInd[i]] .= i
    end

    idx = copy(idxNew)

    if ( totalProcesses > 1)
        idxCount = countmap(idx)
        countVector = [get(idxCount, i, 0) for i in 1:maximum(idx)]
        processStab = silhouettes(trunc.(Int, idx), countVector, pairwise(CosineDist(), Wall, dims=2))
        processStabAvg = zeros(1, totalProcesses)
        for i = 1 : totalProcesses
            processStabAvg[i] = mean(processStab[idx.==i])
        end
    else
        allDist = pairwise(CosineDist(), hcat(centroids', Wall), dims=2)
        processStab = 1 - allDist[1:size(centroids', 2), (size(centroids', 2)+1): size(allDist, 2) ]'
        processStabAvg = mean(processStab)
    end

    centroidStd = zeros( size(centroids) )
    for i = 1 : totalProcesses
        centroidStd[i,:] = std(Wall[:, idx.==i]; dims=2 )
    end

    centroids = centroids'
    centroidStd = centroidStd'

    idxS = zeros(size(idx))
    for i = 1 : totalProcesses : size(Wall, 2)
        iEnd = i + totalProcesses - 1
        idxG = idx[i:iEnd]

        for j = 1 : totalProcesses
            idxS[i+j-1] = findfirst(idxG .== j)
        end
    end

    exposure = zeros(convert(Int, maximum(idxS)), size(Hall, 2))
    exposureStd = zeros(size(exposure))

    for i = 1 : convert(Int, maximum(idxS))
        exposure[i, :] = mean(Hall[idx.==i, :], dims=1)
        exposureStd[i, :] = std(Hall[idx.==i, :], dims=1)
    end

    return centroids, centroidStd, exposure, exposureStd, idx, idxS, processStab, processStabAvg, clusterCompactness
end

max_stars_repo_name: oriolpich/sigprofilerjulia

max_issues_repo_name: oriolpich/sigprofilerjulia

max_forks_repo_name: oriolpich/sigprofilerjulia


text:
using CalculusWithJulia
using Test

include("package-test.jl")
!isinteractive() && include("doc-build-test.jl")


text:
include("EllCrv/EllCrv.jl")
include("EllCrv/Finite.jl")
include("EllCrv/Heights.jl")
include("EllCrv/MinimalModels.jl")
include("EllCrv/Misc.jl")
include("EllCrv/Models.jl")
include("EllCrv/MordellWeilRank.jl")
include("EllCrv/Isogeny.jl")
include("EllCrv/Isomorphisms.jl")
include("EllCrv/Torsion.jl")
include("EllCrv/Twists.jl")


text:
## This file contains adapter code for doing depthwise convolutions with im2col.

"""
    depthwiseconv_im2col!(y, x, w, cdims, col=similar(x); alpha=1, beta=0)

Perform a depthwise convolution using im2col and GEMM, store the result in `y`.

See `conv_im2col!()` for an explanation of optional parameters.
"""
depthwiseconv_im2col!

function depthwiseconv_im2col!(
                y::AbstractArray{T,5}, x::AbstractArray{T,5},
                w::AbstractArray{T,5}, cdims::DepthwiseConvDims;
                col::AbstractArray{T,3} = similar(x, im2col_dims(cdims)),
                alpha::T=T(1), beta::T=T(0)) where T
    check_dims(size(x), size(w), size(y), cdims)

    # This functions exactly the same as conv_im2col!(), except that we shard the
    # incoming data into slices of single channels.  This means that we need to walk
    # each pointer forward individually, as done below, taking a single input channel
    # and combining it with each kernel individually, before walking forward and doing
    # the next input channel.
    M = prod(output_size(cdims))
    N = channel_multiplier(cdims)
    K = prod(kernel_size(cdims))

    dcdims = DenseConvDims(cdims)
    @threads for batch_idx in 1:size(x)[end]
        # col_slice is a thread-local workspace
        col_slice = view(col, :, :, threadid())

        im2col!(col_slice, view(x, :, :, :, :, batch_idx), dcdims)

        # We do a separate convolution for each channel in x, as we must
        for c_in in 1:channels_in(cdims)
            # Walk each pointer forward as we process each input channel
            GC.@preserve col_slice w y begin
                col_ptr = pointer(col_slice, (c_in-1)*M*K+1)
                w_ptr = pointer(w, (c_in-1)*K*N+1)
                y_ptr = pointer(y, ((batch_idx - 1)*channels_in(cdims) + c_in - 1)*M*N + 1)
                gemm!(Val(false), Val(false), M, N, K, alpha, col_ptr, w_ptr, beta, y_ptr)
            end
        end
    end
    return y
end

"""
    ∇depthwiseconv_filter_im2col!(dw, w, dy, cdims, col=similar(dw); alpha=1, beta)

Depthwise conv2d backward pass onto the weights using im2col and GEMM.
See the documentation for `conv_im2col!()` for explanation of optional parameters.
"""
∇depthwiseconv_filter_im2col!

function ∇depthwiseconv_filter_im2col!(
                dw::AbstractArray{T,5}, x::AbstractArray{T,5},
                dy::AbstractArray{T,5}, cdims::DepthwiseConvDims;
                col::AbstractArray{T,3} = similar(dw, im2col_dims(cdims)),
                alpha::T=T(1), beta::T=T(0)) where T
    check_dims(size(x), size(dw), size(dy), cdims)

    M = prod(kernel_size(cdims))
    N = channel_multiplier(cdims)
    K = prod(output_size(cdims))

    @threads for batch_idx in 1:size(x)[end]
        # col_slice is a thread-local workspace
        col_slice = view(col, :, :, threadid())
        im2col!(col_slice, view(x, :, :, :, :, batch_idx), cdims)

        # We do a separate convolution for each channel in x, as we must
        for c_in in 1:channels_in(cdims)
            # Walk each pointer forward as we process each input channel
            GC.@preserve col_slice dw dy begin
                col_ptr = pointer(col_slice, (c_in - 1)*M*K + 1)
                dy_ptr = pointer(dy, (batch_idx - 1)*N*K*channels_in(cdims) + (c_in - 1)*K*N + 1)
                dw_ptr = pointer(dw, (c_in - 1)*M*N + 1)
                gemm!(Val(true), Val(false), M, N, K, alpha, col_ptr, dy_ptr, beta, dw_ptr)
            end
        end

        # Because we accumulate over batches in this loop, we must set `beta` equal
        # to `1.0` from this point on.
        beta = T(1)
    end
    return dw
end

"""
    depthwiseconv2d_Δx_im2col!(dx, w, dy, cdims, col=similar(dx); alpha=1, beta=0)

Depwthwise conv2d backward pass onto the input using im2col and GEMM.
See the documentation for `conv_im2col!()` for explanation of optional parameters.
"""
∇depthwiseconv_data_im2col!

function ∇depthwiseconv_data_im2col!(
                dx::AbstractArray{T,5}, dy::AbstractArray{T,5},
                w::AbstractArray{T,5}, cdims::DepthwiseConvDims;
                col::AbstractArray{T,3} = similar(dx, im2col_dims(cdims)),
                alpha::T=T(1), beta::T=T(0)) where T
    check_dims(size(dx), size(w), size(dy), cdims)

    M = prod(output_size(cdims))
    N = prod(kernel_size(cdims))
    K = channel_multiplier(cdims)

    @threads for batch_idx in 1:size(dx)[end]
        # col_slice is a thread-local workspace
        col_slice = view(col, :, :, threadid())

        # We do a separate convolution for each channel in x, as we must
        for cidx in 1:channels_in(cdims)
            GC.@preserve col_slice w dy begin
                # Walk each pointer forward as we process each input channel
                dy_ptr = pointer(dy, (batch_idx - 1)*M*K*channels_in(cdims)+(cidx - 1)*K*M + 1)
                w_ptr = pointer(w, (cidx - 1)*K*N + 1)
                col_ptr = pointer(col_slice, (cidx - 1)*M*N + 1)
                gemm!(Val(false), Val(true), M, N, K, alpha, dy_ptr, w_ptr, T(0), col_ptr)
            end
        end
        col2im!(view(dx, :, :, :, :, batch_idx), col_slice, cdims)
    end
    return dx
end


text:
module alsa
	module libasound
		using sys

		@pkgconf begin
			using libc.POSIX

			c`-lasound`
			c"""
				#include <alsa/asoundlib.h>
			"""ji
		end
	end
end


text:
@symbol_func function P_LH(cur_reactor::AbstractReactor)
  cur_P = 0.0976

  cur_P *= cur_reactor.n_bar ^ 0.72

  cur_P *= cur_reactor.B_0 ^ 0.80

  cur_P *= surface_area(cur_reactor) ^ 0.94

  cur_P /= cur_reactor.A

  cur_P
end


text:
# ---
# title: 1438. Longest Continuous Subarray With Absolute Diff Less Than or Equal to Limit
# id: problem1438
# author: Qling
# date: 2021-02-21
# difficulty: Medium
# categories: Array, Sliding Window
# link: <https://leetcode.com/problems/longest-continuous-subarray-with-absolute-diff-less-than-or-equal-to-limit/description/>
# hidden: true
# ---
#
# Given an array of integers `nums` and an integer `limit`, return the size of
# the longest **non-empty** subarray such that the absolute difference between
# any two elements of this subarray is less than or equal to `limit` _._
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [8,2,4,7], limit = 4
#     Output: 2
#     Explanation: All subarrays are:
#     [8] with maximum absolute diff |8-8| = 0 <= 4.
#     [8,2] with maximum absolute diff |8-2| = 6 > 4.
#     [8,2,4] with maximum absolute diff |8-2| = 6 > 4.
#     [8,2,4,7] with maximum absolute diff |8-2| = 6 > 4.
#     [2] with maximum absolute diff |2-2| = 0 <= 4.
#     [2,4] with maximum absolute diff |2-4| = 2 <= 4.
#     [2,4,7] with maximum absolute diff |2-7| = 5 > 4.
#     [4] with maximum absolute diff |4-4| = 0 <= 4.
#     [4,7] with maximum absolute diff |4-7| = 3 <= 4.
#     [7] with maximum absolute diff |7-7| = 0 <= 4.
#     Therefore, the size of the longest subarray is 2.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [10,1,2,4,7,2], limit = 5
#     Output: 4
#     Explanation: The subarray [2,4,7,2] is the longest since the maximum absolute diff is |2-7| = 5 <= 5.
#
#
# **Example 3:**
#
#
#
#     Input: nums = [4,2,2,2,4,4,2,2], limit = 0
#     Output: 3
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 10^5`
#   * `1 <= nums[i] <= 10^9`
#   * `0 <= limit <= 10^9`
#
#
## @lc code=start

using LeetCode

function longest_subarray(nums::Vector{Int}, limit::Int)
    queMin, queMax = Deque{Int}(), Deque{Int}()
    n = length(nums)
    left = right = 1
    res = 0

    while right <= n
        while !isempty(queMin) && !isempty(queMax) && last(queMax) < nums[right]
            pop!(queMax)
        end
        while !isempty(queMin) && !isempty(queMax) && last(queMin) > nums[right]
            pop!(queMin)
        end
        push!(queMax, nums[right])
        push!(queMin, nums[right])

        while !isempty(queMin) && !isempty(queMax) && first(queMax) - first(queMin) > limit
            if nums[left] == first(queMin)
                popfirst!(queMin)
            end
            if nums[left] == first(queMax)
                popfirst!(queMax)
            end
            left += 1
        end

        res = max(res, right - left + 1)
        right += 1

    end

    return res
end
## @lc code=end


text:
include("../emi.jl")
using EMI
using EMI.Draw, EMI.Gmsh

const R = 1
const C = R/sqrt(2)

loop = Loop([CircleArc(Point(-R, 0), Point(0, 0), Point(0, R)),
             CircleArc(Point(0, R), Point(0, 0), Point(C, C)),
             Line(Point(C, C), Point(0., 0.)),
             Line(Point(0., 0.), Point(C, -C)),
             CircleArc(Point(C, -C), Point(0., 0.), Point(0, -R)),
             CircleArc(Point(0, -R), Point(0, 0), Point(-R, 0))])

canvas = Canvas()
canvas = canvas + loop

set_bbox!(canvas, 0.2, 0.2)
gmsh_script(canvas, 0.2)

max_stars_repo_name: MiroK/emi-cylinders

max_issues_repo_name: MiroK/emi-cylinders

max_forks_repo_name: MiroK/emi-cylinders


text:
mutable struct ReadCoverageAccumulator
	dCoverage::Vector{Vector{Int}} # derivative of coverage - makes accumulation faster since only starting/ending points of read chunks must be collected
end

function ReadCoverageAccumulator(seqLengths::Vector{Int})
	ReadCoverageAccumulator(Vector{Int}[zeros(Int,len+1) for len in seqLengths]) # length + 1 to make sure we end at 0
end

function _integratecoverage(x)
	y = cumsum(x)
	@assert y[end]==0
	y[1:end-1]
end

function coverage(rca::ReadCoverageAccumulator)
	N = length(rca.dCoverage)
	map!(_integratecoverage, Vector{Vector{Int}}(undef,N), rca.dCoverage)
end

function processread!(rca::ReadCoverageAccumulator, r::BamRead)
	refInd = refID(r)+1 # 0 to 1-based index
	dCoverage = rca.dCoverage[refInd]

	for c in r
		o = op(c)

		if o==CIGAR_M || o==CIGAR_X || o==CIGAR_Eq # the only chunks that are aligned to both read and reference
			pStart = refpos(c)
			dCoverage[pStart] += 1

			pEnd = pStart+len(c)
			dCoverage[pEnd] -= 1
		end
	end
end

function readcoverage(bamFile::BamFile; strands=:both, mappingQualityThreshold::Int=30, ignoreChimeric::Bool=true)
	@assert strands∈[:both, :forward, :reverse] "Strand must be one of :both, :forward, :reverse."

	seqs = sequences(bamFile)
	seqLengths = Int[s[2] for s in seqs]

	flagsUnset = BAM_FUNMAP | BAM_FSECONDARY | BAM_FDUP
	flagsSet = 0

	strands == :forward && (flagsUnset |= BAM_FREVERSE)
	strands == :reverse && (flagsSet   |= BAM_FREVERSE)

	rca = ReadCoverageAccumulator(seqLengths)
	processbam!(bamFile, rca, flagsSet=flagsSet, flagsUnset=flagsUnset,
		        mappingQualityThreshold=mappingQualityThreshold,
		        ignoreChimeric=ignoreChimeric)

	coverage(rca)
end

readcoverage(bamFilename::AbstractString; kwargs...) = readcoverage(BamFile(bamFilename);kwargs...)

# get a read coverage matrix per reference sequence
function readcoverage(inputs::AbstractArray; kwargs...)
	nbrSamples = length(inputs)

	coverage = map(x->readcoverage(x;kwargs...),inputs) # Vector{Vector{Vector{Int}}} : Samples-ReferenceSequence-Coverage

	# reorganize to create one matrix per reference sequence
	nbrSeqsPerSample = map(length, coverage)
	nbrSeqs = nbrSeqsPerSample[1]

	@assert all(nbrSeqsPerSample .== nbrSeqs) "All samples must have the same number of reference sequences."

	covMatrices = Vector{Matrix{Int}}(undef,nbrSeqs)
	for i=1:nbrSeqs
		seqLengthPerSample = map( x->length(x[i]), coverage )
		seqLength = maximum(seqLengthPerSample) # To avoid failing if some samples have indels

		# gather read coverage in matrix filled with zeros at the end
		M = zeros(nbrSamples, seqLength)
		for j=1:nbrSamples
			s = seqLengthPerSample[j]
			M[j,1:s] = coverage[j][i]
		end
		covMatrices[i] = M
	end
	covMatrices
end


text:
""" Solve the classical GP Regression (especially valid for augmented likelihoods) """
mutable struct Analytic{T<:Real} <: Inference{T}
    ϵ::T #Convergence criteria
    nIter::Integer #Number of steps performed
    Stochastic::Bool #Use of mini-batches
    nSamples::Int64
    nSamplesUsed::Int64 #Size of mini-batches
    MBIndices::Vector{Int64} #Indices of the minibatch
    ρ::T #Stochastic Coefficient
    HyperParametersUpdated::Bool #To know if the inverse kernel matrix must updated
    function Analytic{T}(ϵ::T,nIter::Integer,Stochastic::Bool,nSamples::Integer,MBIndices::AbstractVector,ρ::T,flag::Bool) where T
        return new{T}(ϵ,nIter,Stochastic,nSamples,nSamples,MBIndices,ρ,flag)
    end
end

"""`Analytic(;ϵ::T=1e-5)`

Analytic inference structure for the classical GP regression

**Keyword arguments**

    - `ϵ::T` : convergence criteria, which can be user defined
"""
function Analytic(;ϵ::T=1e-5) where {T<:Real}
    Analytic{T}(ϵ,0,false,1,collect(1:1),1.0,true)
end

function Base.show(io::IO,inference::Analytic{T}) where T
    print(io,"Analytic Inference")
end

function init_inference(inference::Analytic{T},nLatent::Integer,nFeatures::Integer,nSamples::Integer,nSamplesUsed::Integer) where {T<:Real}
    inference.nSamples = nSamples
    inference.nSamplesUsed = nSamples
    inference.MBIndices = 1:nSamples
    inference.ρ = nSamples/nSamplesUsed
    return inference
end


text:
# to test vtk-files
using StableRNGs
OVERWRITE_CHECKSUMS = false
checksums_file = joinpath(dirname(@__FILE__), "checksums.sha1")
checksum_list = read(checksums_file, String)
if OVERWRITE_CHECKSUMS
    csio = open(checksums_file, "w")
else
    csio = open(checksums_file, "r")
end

@testset "Grid, DofHandler, vtk" begin
    for (celltype, dim) in ((Line,                   1),
                            (QuadraticLine,          1),
                            (Quadrilateral,          2),
                            (QuadraticQuadrilateral, 2),
                            (Triangle,               2),
                            (QuadraticTriangle,      2),
                            (Hexahedron,             3),
                            (Cell{3,20,6},           3),
                            (Tetrahedron,            3))

        # create test grid, do some operations on it and then test
        # the resulting sha1 of the stored vtk file
        # after manually checking the exported vtk
        nels = ntuple(x->5, dim)
        right = Vec{dim, Float64}(ntuple(x->1.5, dim))
        left = -right
        grid = generate_grid(celltype, nels, left, right)

        transform!(grid, x-> 2x)

        radius = 2*1.5
        addcellset!(grid, "cell-1", [1,])
        addcellset!(grid, "middle-cells", x -> norm(x) < radius)
        addfaceset!(grid, "middle-faceset", x -> norm(x) < radius)
        addfaceset!(grid, "right-faceset", getfaceset(grid, "right"))
        addnodeset!(grid, "middle-nodes", x -> norm(x) < radius)

        gridfilename = "grid-$(Ferrite.celltypes[celltype])"
        vtk_grid(gridfilename, grid) do vtk
            vtk_cellset(vtk, grid, "cell-1")
            vtk_cellset(vtk, grid, "middle-cells")
            vtk_nodeset(vtk, grid, "middle-nodes")
        end

        # test the sha of the file
        sha = bytes2hex(open(SHA.sha1, gridfilename*".vtu"))
        if OVERWRITE_CHECKSUMS
            write(csio, sha, "\n")
        else
            @test chomp(readline(csio)) == sha
            rm(gridfilename*".vtu")
        end

        # Create a DofHandler, add some things, write to file and
        # then check the resulting sha
        dofhandler = DofHandler(grid)
        push!(dofhandler, :temperature, 1)
        push!(dofhandler, :displacement, 3)
        close!(dofhandler)
        ch = ConstraintHandler(dofhandler)
        dbc = Dirichlet(:temperature, union(getfaceset(grid, "left"), getfaceset(grid, "right-faceset")), (x,t)->1)
        add!(ch, dbc)
        dbc = Dirichlet(:temperature, getfaceset(grid, "middle-faceset"), (x,t)->4)
        add!(ch, dbc)
        for d in 1:dim
            dbc = Dirichlet(:displacement, union(getfaceset(grid, "left")), (x,t) -> d, d)
            add!(ch, dbc)
        end
        close!(ch)
        update!(ch, 0.0)
        rng = StableRNG(1234)
        u = rand(rng, ndofs(dofhandler))
        apply!(u, ch)

        dofhandlerfilename = "dofhandler-$(Ferrite.celltypes[celltype])"
        vtk_grid(dofhandlerfilename, dofhandler) do vtk
            vtk_point_data(vtk, ch)
            vtk_point_data(vtk, dofhandler, u)
        end

        # test the sha of the file
        sha = bytes2hex(open(SHA.sha1, dofhandlerfilename*".vtu"))
        if OVERWRITE_CHECKSUMS
            write(csio, sha, "\n")
        else
            @test chomp(readline(csio)) == sha
            rm(dofhandlerfilename*".vtu")
        end

    end

end # of testset

close(csio)

@testset "vtk tensor export" begin
    # open files
    checksums_file_tensors = joinpath(dirname(@__FILE__), "checksums2.sha1")
    if OVERWRITE_CHECKSUMS
        csio = open(checksums_file_tensors, "w")
    else
        csio = open(checksums_file_tensors, "r")
    end

    # 3D grid
    grid = generate_grid(Hexahedron, (1,1,1))

    sym_tensor_data = [SymmetricTensor{2,3}(ntuple(i->i, 6)) for j=1.0:8.0]
    tensor_data = [Tensor{2,3}(ntuple(i->i, 9)) for j=1.0:8.0]
    vector_data = [Vec{3}(ntuple(i->i, 3)) for j=1:8]

    filename_3d = "test_vtk_3d"
    vtk_grid(filename_3d, grid) do vtk_file
        vtk_point_data(vtk_file, sym_tensor_data, "symmetric tensor")
        vtk_point_data(vtk_file, tensor_data, "tensor")
        vtk_point_data(vtk_file, vector_data, "vector")
    end

    # 2D grid
    grid = generate_grid(Quadrilateral, (1,1))

    sym_tensor_data = [SymmetricTensor{2,2}(ntuple(i->i, 3)) for j=1.0:4.0]
    tensor_data = [Tensor{2,2}(ntuple(i->i, 4)) for j=1.0:4.0]
    tensor_data_1D = [SymmetricTensor{2,1}(ntuple(i->i, 1)) for j=1.0:4.0]
    vector_data = [Vec{2}(ntuple(i->i, 2)) for j=1:4]

    filename_2d = "test_vtk_2d"
    vtk_grid(filename_2d, grid) do vtk_file
        vtk_point_data(vtk_file, sym_tensor_data, "symmetric tensor")
        vtk_point_data(vtk_file, tensor_data, "tensor")
        vtk_point_data(vtk_file, tensor_data_1D, "tensor_1d")
        vtk_point_data(vtk_file, vector_data, "vector")
    end

    # test the shas of the files
    files = [filename_3d, filename_2d]
    for filename in files
        sha = bytes2hex(open(SHA.sha1, filename*".vtu"))
        if OVERWRITE_CHECKSUMS
            write(csio, sha, "\n")
        else
            @test chomp(readline(csio)) == sha
            rm(filename*".vtu")
        end
    end

    close(csio)
end

@testset "Grid utils" begin

    grid = Ferrite.generate_grid(QuadraticQuadrilateral, (1, 1), Vec((0.,0.)), Vec((1.,1.)))

    addcellset!(grid, "cell_set", [1]);
    node_set = Set(1:getnnodes(grid))
    addnodeset!(grid, "node_set", node_set)

    @test getnodesets(grid) == Dict("node_set" => node_set)

    @test getnodes(grid, [1]) == [getnodes(grid, 1)] # untested

    @test length(getnodes(grid, "node_set")) == 9

    @test collect(getcoordinates(getnodes(grid, 5)).data) ≈ [0.5, 0.5]

    @test getcells(grid, "cell_set") == [getcells(grid, 1)]

    f(x) = Tensor{1,1,Float64}((1 + x[1]^2 + 2x[2]^2, ))

    values = compute_vertex_values(grid, f)
    @test f([0.0, 0.0]) == values[1]
    @test f([0.5, 0.5]) == values[5]
    @test f([1.0, 1.0]) == values[9]

    @test compute_vertex_values(grid, collect(1:9), f) == values

    # Can we test this in a better way? The set makes the order random.
    @test length(compute_vertex_values(grid, "node_set", f)) == 9

    # CellIterator on a grid without DofHandler
    grid = generate_grid(Triangle, (4,4))
    n = 0
    ci = CellIterator(grid)
    @test length(ci) == getncells(grid)
    for c in ci
        getcoordinates(c)
        getnodes(c)
        n += cellid(c)
    end
    @test n == div(getncells(grid)*(getncells(grid) + 1), 2)
end

@testset "Grid sets" begin

    grid = Ferrite.generate_grid(Hexahedron, (1, 1, 1), Vec((0.,0., 0.)), Vec((1.,1.,1.)))

    #Test manual add
    addcellset!(grid, "cell_set", [1]);
    addnodeset!(grid, "node_set", [1])
    addfaceset!(grid, "face_set", [FaceIndex(1,1)])
    addedgeset!(grid, "edge_set", [EdgeIndex(1,1)])
    addvertexset!(grid, "vert_set", [VertexIndex(1,1)])

    #Test function add
    addfaceset!(grid, "left_face", (x)-> x[1] ≈ 0.0)
    addedgeset!(grid, "left_lower_edge", (x)-> x[1] ≈ 0.0 && x[3] ≈ 0.0)
    addvertexset!(grid, "left_corner", (x)-> x[1] ≈ 0.0 && x[2] ≈ 0.0 && x[3] ≈ 0.0)

    @test 1 in Ferrite.getnodeset(grid, "node_set")
    @test FaceIndex(1,5) in getfaceset(grid, "left_face")
    @test EdgeIndex(1,4) in getedgeset(grid, "left_lower_edge")
    @test VertexIndex(1,1) in getvertexset(grid, "left_corner")

end

@testset "Grid topology" begin
#                           (11)
#                   (10)+-----+-----+(12)
#                       |  5  |  6  |
#                   (7) +-----+-----+(9)
#                       |  3  |  4  |
#                   (4) +-----+-----+(6)
#                       |  1  |  2  |
#                   (1) +-----+-----+(3)
#                            (2)
    quadgrid = generate_grid(Quadrilateral,(2,3))
    topology = ExclusiveTopology(quadgrid)
    #test vertex neighbors maps cellid and local vertex id to neighbor id and neighbor local vertex id
    @test topology.vertex_neighbor[1,3] == Ferrite.EntityNeighborhood(VertexIndex(4,1))
    @test topology.vertex_neighbor[2,4] == Ferrite.EntityNeighborhood(VertexIndex(3,2))
    @test topology.vertex_neighbor[3,3] == Ferrite.EntityNeighborhood(VertexIndex(6,1))
    @test topology.vertex_neighbor[3,2] == Ferrite.EntityNeighborhood(VertexIndex(2,4))
    @test topology.vertex_neighbor[4,1] == Ferrite.EntityNeighborhood(VertexIndex(1,3))
    @test topology.vertex_neighbor[4,4] == Ferrite.EntityNeighborhood(VertexIndex(5,2))
    @test topology.vertex_neighbor[5,2] == Ferrite.EntityNeighborhood(VertexIndex(4,4))
    @test topology.vertex_neighbor[6,1] == Ferrite.EntityNeighborhood(VertexIndex(3,3))
    #test face neighbor maps cellid and local face id to neighbor id and neighbor local face id
    @test topology.face_neighbor[1,2] == Ferrite.EntityNeighborhood(FaceIndex(2,4))
    @test topology.face_neighbor[1,3] == Ferrite.EntityNeighborhood(FaceIndex(3,1))
    @test topology.face_neighbor[2,3] == Ferrite.EntityNeighborhood(FaceIndex(4,1))
    @test topology.face_neighbor[2,4] == Ferrite.EntityNeighborhood(FaceIndex(1,2))
    @test topology.face_neighbor[3,1] == Ferrite.EntityNeighborhood(FaceIndex(1,3))
    @test topology.face_neighbor[3,2] == Ferrite.EntityNeighborhood(FaceIndex(4,4))
    @test topology.face_neighbor[3,3] == Ferrite.EntityNeighborhood(FaceIndex(5,1))
    @test topology.face_neighbor[4,1] == Ferrite.EntityNeighborhood(FaceIndex(2,3))
    @test topology.face_neighbor[4,3] == Ferrite.EntityNeighborhood(FaceIndex(6,1))
    @test topology.face_neighbor[4,4] == Ferrite.EntityNeighborhood(FaceIndex(3,2))
    @test topology.face_neighbor[5,1] == Ferrite.EntityNeighborhood(FaceIndex(3,3))
    @test topology.face_neighbor[5,2] == Ferrite.EntityNeighborhood(FaceIndex(6,4))
    @test topology.face_neighbor[5,3] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[])
    @test topology.face_neighbor[5,4] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[])
    @test topology.face_neighbor[6,1] == Ferrite.EntityNeighborhood(FaceIndex(4,3))
    @test topology.face_neighbor[6,2] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[])
    @test topology.face_neighbor[6,3] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[])
    @test topology.face_neighbor[6,4] == Ferrite.EntityNeighborhood(FaceIndex(5,2))
#                         (8)
#                (7) +-----+-----+(9)
#                    |  3  |  4  |
#                (4) +-----+-----+(6) bottom view
#                    |  1  |  2  |
#                (1) +-----+-----+(3)
#                         (2)
#                         (15)
#               (16) +-----+-----+(17)
#                    |  3  |  4  |
#               (13) +-----+-----+(15) top view
#                    |  1  |  2  |
#               (10) +-----+-----+(12)
#                        (11)
    hexgrid = generate_grid(Hexahedron,(2,2,1))
    topology = ExclusiveTopology(hexgrid)
    @test topology.edge_neighbor[1,11] == Ferrite.EntityNeighborhood(EdgeIndex(4,9))
    @test topology.edge_neighbor[2,12] == Ferrite.EntityNeighborhood(EdgeIndex(3,10))
    @test topology.edge_neighbor[3,10] == Ferrite.EntityNeighborhood(EdgeIndex(2,12))
    @test topology.edge_neighbor[4,9] == Ferrite.EntityNeighborhood(EdgeIndex(1,11))
    @test all(iszero,topology.vertex_neighbor)
    @test topology.face_neighbor[1,3] == Ferrite.EntityNeighborhood(FaceIndex(2,5))
    @test topology.face_neighbor[1,4] == Ferrite.EntityNeighborhood(FaceIndex(3,2))
    @test topology.face_neighbor[2,4] == Ferrite.EntityNeighborhood(FaceIndex(4,2))
    @test topology.face_neighbor[2,5] == Ferrite.EntityNeighborhood(FaceIndex(1,3))
    @test topology.face_neighbor[3,2] == Ferrite.EntityNeighborhood(FaceIndex(1,4))
    @test topology.face_neighbor[3,3] == Ferrite.EntityNeighborhood(FaceIndex(4,5))
    @test topology.face_neighbor[4,2] == Ferrite.EntityNeighborhood(FaceIndex(2,4))
    @test topology.face_neighbor[4,5] == Ferrite.EntityNeighborhood(FaceIndex(3,3))

#                   +-----+-----+
#                   |\  6 |\  8 |
#                   |  \  |  \  |
#                   |  5 \| 7  \|
#                   +-----+-----+
#                   |\  2 |\  4 |
#                   |  \  |  \  |
#                   |  1 \| 3  \|
#                   +-----+-----+
# test for multiple vertex_neighbors as in e.g. ele 3, local vertex 3 (middle node)
    trigrid = generate_grid(Triangle,(2,2))
    topology = ExclusiveTopology(trigrid)
    @test topology.vertex_neighbor[3,3] == Ferrite.EntityNeighborhood([VertexIndex(5,2),VertexIndex(6,1),VertexIndex(7,1)])

# test mixed grid
    cells = [
        Hexahedron((1, 2, 3, 4, 5, 6, 7, 8)),
        Quadrilateral((3, 2, 9, 10)),
        ]
    nodes = [Node(coord) for coord in zeros(Vec{2,Float64}, 10)]
    grid = Grid(cells, nodes)
    topology = ExclusiveTopology(grid)
    @test all(iszero,topology.vertex_neighbor)
    @test topology.face_neighbor[2,1] == Ferrite.EntityNeighborhood(EdgeIndex(1,2))
    @test topology.edge_neighbor[1,2] == Ferrite.EntityNeighborhood(FaceIndex(2,1))
#
#                   +-----+-----+-----+
#                   |  7  |  8  |  9  |
#                   +-----+-----+-----+
#                   |  4  |  5  |  6  |
#                   +-----+-----+-----+
#                   |  1  |  2  |  3  |
#                   +-----+-----+-----+
# test application: form level 1 neighborhood patches of elements
    quadgrid = generate_grid(Quadrilateral,(3,3))
    topology = ExclusiveTopology(quadgrid)
    patches = Vector{Int}[Ferrite.getneighborhood(topology, quadgrid, CellIndex(i)) for i in 1:getncells(quadgrid)]

    @test issubset([4,5,2], patches[1]) # neighbor elements of element 1 are 4 5 and 2
    @test issubset([1,4,5,6,3], patches[2])
    @test issubset([2,5,6], patches[3])
    @test issubset([7,8,5,2,1], patches[4])
    @test issubset([1,2,3,4,6,7,8,9], patches[5])
    @test issubset([3,2,5,8,9], patches[6])
    @test issubset([4,5,8], patches[7])
    @test issubset([7,4,5,6,9], patches[8])
    @test issubset([8,5,6], patches[9])

    @test Ferrite.getneighborhood(topology, quadgrid, VertexIndex(1,1)) == [VertexIndex(1,2), VertexIndex(1,4)]
    @test Ferrite.getneighborhood(topology, quadgrid, VertexIndex(2,1)) == [VertexIndex(1,1), VertexIndex(1,3), VertexIndex(2,2), VertexIndex(2,4)]
    @test Ferrite.getneighborhood(topology, quadgrid, VertexIndex(5,4)) == [VertexIndex(4,2), VertexIndex(4,4), VertexIndex(5,1), VertexIndex(5,3), VertexIndex(7,1), VertexIndex(7,3), VertexIndex(8,2), VertexIndex(8,4)]
    @test Ferrite.toglobal(quadgrid, Ferrite.getneighborhood(topology, quadgrid, VertexIndex(1,1))) == [2,5]
    @test Ferrite.toglobal(quadgrid, Ferrite.getneighborhood(topology, quadgrid, VertexIndex(2,1))) == [1,6,3]
    @test Ferrite.toglobal(quadgrid, Ferrite.getneighborhood(topology, quadgrid, VertexIndex(5,4))) == [6,9,11,14]

    @test topology.face_skeleton == [FaceIndex(1,1),FaceIndex(1,2),FaceIndex(1,3),FaceIndex(1,4),
                                          FaceIndex(2,1),FaceIndex(2,2),FaceIndex(2,3),
                                          FaceIndex(3,1),FaceIndex(3,2),FaceIndex(3,3),
                                          FaceIndex(4,2),FaceIndex(4,3),FaceIndex(4,4),
                                          FaceIndex(5,2),FaceIndex(5,3),FaceIndex(6,2),FaceIndex(6,3),
                                          FaceIndex(7,2),FaceIndex(7,3),FaceIndex(7,4),
                                          FaceIndex(8,2),FaceIndex(8,3),FaceIndex(9,2),FaceIndex(9,3)]
    @test length(topology.face_skeleton) == 4*3 + 3*4

    quadratic_quadgrid = generate_grid(QuadraticQuadrilateral,(3,3))
    quadgrid_topology = ExclusiveTopology(grid)
    quadgrid_topology.face_skeleton == topology.face_skeleton
#
#                   +-----+-----+-----+
#                   |  7  |  8  |  9  |
#                   +-----+-----+-----+
#                   |  4  |  5  |  6  |
#                   +-----+-----+-----+
#                   |  1  |  2  |  3  |
#                   +-----+-----+-----+
# test application: integrate jump across element boundary 5
    function reinit!(fv::FaceValues, cellid::Int, faceid::Int, grid)
        coords = getcoordinates(grid, cellid)
        Ferrite.reinit!(fv, coords, faceid)
    end
    reinit!(fv::FaceValues, faceid::FaceIndex, grid) = reinit!(fv,faceid[1],faceid[2],grid) # wrapper for reinit!(fv,cellid,faceid,grid)
    face_neighbors_ele5 = nonzeros(topology.face_neighbor[5,:])
    ip = Lagrange{2, RefCube, 1}()
    qr_face = QuadratureRule{1, RefCube}(2)
    fv_ele = FaceVectorValues(qr_face, ip)
    fv_neighbor = FaceVectorValues(qr_face, ip)
    u_ele5 = [3.0 for _ in 1:8]
    u_neighbors = [5.0 for _ in 1:8]
    jump_int = 0.
    jump_abs = 0.
    for ele_faceid in 1:nfaces(quadgrid.cells[5])
        reinit!(fv_ele, 5, ele_faceid, quadgrid)
        for q_point in 1:getnquadpoints(fv_ele)
            dΩ = getdetJdV(fv_ele, q_point)
            normal_5 = getnormal(fv_ele, q_point)
            u_5_n = function_value(fv_ele, q_point, u_ele5) ⋅ normal_5
            for neighbor_entity in face_neighbors_ele5[ele_faceid].neighbor_info # only one entity can be changed to get rid of the for loop
                reinit!(fv_neighbor, neighbor_entity, quadgrid)
                normal_neighbor = getnormal(fv_neighbor, q_point)
                u_neighbor = function_value(fv_neighbor, q_point, u_neighbors) ⋅ normal_neighbor
                jump_int += (u_5_n + u_neighbor) * dΩ
                jump_abs += abs(u_5_n + u_neighbor) * dΩ
            end
        end
    end
    @test isapprox(jump_abs, 2/3*2*4,atol=1e-6) # 2*4*0.66666, jump is always 2, 4 sides, length =0.66
    @test isapprox(jump_int, 0.0, atol=1e-6)
end

@testset "grid coloring" begin
    function test_coloring(grid, cellset=Set(1:getncells(grid)))
        for alg in (ColoringAlgorithm.Greedy, ColoringAlgorithm.WorkStream)
            color_vectors = create_coloring(grid, cellset; alg=alg)
            @test sum(length, color_vectors) == length(cellset)
            @test union(Set.(color_vectors)...) == cellset
            conn = Ferrite.create_incidence_matrix(grid, cellset)
            for color in color_vectors, c1 in color, c2 in color
                @test !conn[c1, c2]
            end
        end
    end
    test_coloring(generate_grid(Line, (5,)))
    test_coloring(generate_grid(QuadraticLine, (5,)))
    test_coloring(generate_grid(Triangle, (5, 5)))
    test_coloring(generate_grid(QuadraticTriangle, (5, 5)))
    test_coloring(generate_grid(Quadrilateral, (5, 5)))
    test_coloring(generate_grid(QuadraticQuadrilateral, (5, 5)))
    test_coloring(generate_grid(Tetrahedron, (5, 5, 5)))
    # test_coloring(generate_grid(QuadraticTetrahedron, (5, 5, 5)))
    test_coloring(generate_grid(Hexahedron, (5, 5, 5)))
    # test_coloring(generate_grid(QuadraticHexahedron, (5, 5, 5)))

    # color only a subset
    test_coloring(generate_grid(Line, (5,)), Set{Int}(1:3))
    test_coloring(generate_grid(Triangle, (5, 5)), Set{Int}(1:3^2))
    test_coloring(generate_grid(Quadrilateral, (5, 5)), Set{Int}(1:3^2))
    test_coloring(generate_grid(Tetrahedron, (5, 5, 5)), Set{Int}(1:3^3))
    test_coloring(generate_grid(Hexahedron, (5, 5, 5)), Set{Int}(1:3^3))
    # unconnected subset
    test_coloring(generate_grid(Triangle, (10, 10)), union(Set(1:10), Set(70:80)))
end


text:
module ArrayInterface

using Requires
using LinearAlgebra
using SparseArrays

function ismutable end

"""
    ismutable(x::DataType)

Query whether a type is mutable or not, see
https://github.com/JuliaDiffEq/RecursiveArrayTools.jl/issues/19.
"""
Base.@pure ismutable(x::DataType) = x.mutable
ismutable(x) = ismutable(typeof(x))

ismutable(::Type{<:Array}) = true
ismutable(::Type{<:Number}) = false

# Piracy
function Base.setindex(x::Array,v,i::Int)
  _x = copy(x)
  _x[i] = v
  _x
end

"""
    can_setindex(x::DataType)

Query whether a type can use `setindex!`
"""
can_setindex(x) = true
can_setindex(x::AbstractArray) = can_setindex(typeof(x))

"""
    fast_scalar_indexing(x)

Query whether an array type has fast scalar indexing
"""
fast_scalar_indexing(x) = true
fast_scalar_indexing(x::AbstractArray) = fast_scalar_indexing(typeof(x))

"""
    allowed_getindex(x,i...)

A scalar getindex which is always allowed
"""
allowed_getindex(x,i...) = x[i...]

"""
    allowed_setindex!(x,v,i...)

A scalar setindex! which is always allowed
"""
allowed_setindex!(x,v,i...) = setindex!(x,v,i...)

"""
    isstructured(x::DataType)

Query whether a type is a representation of a structured matrix
"""
isstructured(x) = false
isstructured(x::AbstractArray) = isstructured(typeof(x))
isstructured(::Symmetric) = true
isstructured(::Hermitian) = true
isstructured(::UpperTriangular) = true
isstructured(::LowerTriangular) = true
isstructured(::Tridiagonal) = true
isstructured(::SymTridiagonal) = true
isstructured(::Bidiagonal) = true
isstructured(::Diagonal) = true

"""
    has_sparsestruct(x::AbstractArray)

determine whether `findstructralnz` accepts the parameter `x`
"""
has_sparsestruct(x) = false
has_sparsestruct(x::AbstractArray) = has_sparsestruct(typeof(x))
has_sparsestruct(x::Type{<:AbstractArray}) = false
has_sparsestruct(x::Type{<:SparseMatrixCSC}) = true
has_sparsestruct(x::Type{<:Diagonal}) = true
has_sparsestruct(x::Type{<:Bidiagonal}) = true
has_sparsestruct(x::Type{<:Tridiagonal}) = true
has_sparsestruct(x::Type{<:SymTridiagonal}) = true

"""
    findstructralnz(x::AbstractArray)

Return: (I,J) #indexable objects
Find sparsity pattern of special matrices, the same as the first two elements of findnz(::SparseMatrixCSC)
"""
function findstructralnz(x::Diagonal)
  n=size(x,1)
  (1:n,1:n)
end

abstract type MatrixIndex end

struct BidiagonalIndex <: MatrixIndex
  count::Int
  isup::Bool
end

struct TridiagonalIndex <: MatrixIndex
  count::Int
  nsize::Int
  isrow::Bool
end

Base.firstindex(ind::MatrixIndex)=1
Base.lastindex(ind::MatrixIndex)=ind.count
Base.length(ind::MatrixIndex)=ind.count
function Base.getindex(ind::BidiagonalIndex,i::Int)
  1 <= i <= ind.count || throw(BoundsError(ind, i))
  if ind.isup
    ii=i+1
  else
    ii=i+1+1
  end
  convert(Int,floor(ii/2))
end

function Base.getindex(ind::TridiagonalIndex,i::Int)
  1 <= i <= ind.count || throw(BoundsError(ind, i))
  offsetu= ind.isrow ? 0 : 1
  offsetl= ind.isrow ? 1 : 0
  if 1 <= i <= ind.nsize
    return i
  elseif ind.nsize < i <= ind.nsize+ind.nsize-1
    return i-ind.nsize+offsetu
  else
    return i-(ind.nsize+ind.nsize-1)+offsetl
  end
end

function findstructralnz(x::Bidiagonal)
  n=size(x,1)
  isup= x.uplo=='U' ? true : false
  rowind=BidiagonalIndex(n+n-1,isup)
  colind=BidiagonalIndex(n+n-1,!isup)
  (rowind,colind)
end

function findstructralnz(x::Union{Tridiagonal,SymTridiagonal})
  n=size(x,1)
  rowind=TridiagonalIndex(n+n-1+n-1,n,true)
  colind=TridiagonalIndex(n+n-1+n-1,n,false)
  (rowind,colind)
end

function findstructralnz(x::SparseMatrixCSC)
  rowind,colind,_=findnz(x)
  (rowind,colind)
end

abstract type ColoringAlgorithm end

"""
    fast_matrix_colors(A)

    Query whether a matrix has a fast algorithm for getting the structural
    colors of the matrix.
"""
fast_matrix_colors(A) = false
fast_matrix_colors(A::AbstractArray) = fast_matrix_colors(typeof(A))
fast_matrix_colors(A::Type{<:Union{Diagonal,Bidiagonal,Tridiagonal,SymTridiagonal}}) = true

"""
    matrix_colors(A::Union{Array,UpperTriangular,LowerTriangular})

    The color vector for dense matrix and triangular matrix is simply
    `[1,2,3,...,size(A,2)]`
"""
function matrix_colors(A::Union{Array,UpperTriangular,LowerTriangular})
    eachindex(1:size(A,2)) # Vector size matches number of rows
end

function _cycle(repetend,len)
    repeat(repetend,div(len,length(repetend))+1)[1:len]
end

function matrix_colors(A::Diagonal)
    fill(1,size(A,2))
end

function matrix_colors(A::Bidiagonal)
    _cycle(1:2,size(A,2))
end

function matrix_colors(A::Union{Tridiagonal,SymTridiagonal})
    _cycle(1:3,size(A,2))
end

function __init__()

  @require StaticArrays="90137ffa-7385-5640-81b9-e52037218182" begin
    ismutable(::Type{<:StaticArrays.StaticArray}) = false
    can_setindex(::Type{<:StaticArrays.StaticArray}) = false
    ismutable(::Type{<:StaticArrays.MArray}) = true
  end

  @require LabelledArrays="2ee39098-c373-598a-b85f-a56591580800" begin
    ismutable(::Type{<:LabelledArrays.LArray{T,N,Syms}}) where {T,N,Syms} = ismutable(T)
    can_setindex(::Type{<:LabelledArrays.SLArray}) = false
  end

  @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
    ismutable(::Type{<:Tracker.TrackedArray}) = false
    can_setindex(::Type{<:Tracker.TrackedArray}) = false
    fast_scalar_indexing(::Type{<:Tracker.TrackedArray}) = false
  end

  @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
    fast_scalar_indexing(::Type{<:CuArrays.CuArray}) = false
    @inline allowed_getindex(x::CuArrays.CuArray,i...) = CuArrays._getindex(x,i...)
    @inline allowed_setindex!(x::CuArrays.CuArray,v,i...) = CuArrays._setindex!(x,v,i...)

    function Base.setindex(x::CuArrays.CuArray,v,i::Int)
      _x = copy(x)
      allowed_setindex!(_x,v,i)
      _x
    end
  end

  @require BandedMatrices="aae01518-5342-5314-be14-df237901396f" begin
    is_structured(::Type{<:BandedMatrices.BandedMatrix}) = true
    fast_matrix_colors(::Type{<:BandedMatrices.BandedMatrix}) = true

    function matrix_colors(A::BandedMatrices.BandedMatrix)
        u,l=bandwidths(A)
        width=u+l+1
        _cycle(1:width,size(A,2))
    end

  end

  @require BlockBandedMatrices="aae01518-5342-5314-be14-df237901396f" begin
    is_structured(::Type{<:BandedMatrices.BlockBandedMatrix}) = true
    is_structured(::Type{<:BandedMatrices.BandedBlockBandedMatrix}) = true
    fast_matrix_colors(::Type{<:BlockBandedMatrices.BlockBandedMatrix}) = true
    fast_matrix_colors(::Type{<:BlockBandedMatrices.BandedBlockBandedMatrix}) = true

    function matrix_colors(A::BlockBandedMatrices.BlockBandedMatrix)
        l,u=blockbandwidths(A)
        blockwidth=l+u+1
        nblock=nblocks(A,2)
        cols=[blocksize(A,(1,i))[2] for i in 1:nblock]
        blockcolors=_cycle(1:blockwidth,nblock)
        #the reserved number of colors of a block is the maximum length of columns of blocks with the same block color
        ncolors=[maximum(cols[i:blockwidth:nblock]) for i in 1:blockwidth]
        endinds=cumsum(ncolors)
        startinds=[endinds[i]-ncolors[i]+1 for i in 1:blockwidth]
        colors=[(startinds[blockcolors[i]]:endinds[blockcolors[i]])[1:cols[i]] for i in 1:nblock]
        vcat(colors...)
    end

    function matrix_colors(A::BlockBandedMatrices.BandedBlockBandedMatrix)
        l,u=blockbandwidths(A)
        lambda,mu=subblockbandwidths(A)
        blockwidth=l+u+1
        subblockwidth=lambda+mu+1
        nblock=nblocks(A,2)
        cols=[blocksize(A,(1,i))[2] for i in 1:nblock]
        blockcolors=_cycle(1:blockwidth,nblock)
        #the reserved number of colors of a block is the min of subblockwidth and the largest length of columns of blocks with the same block color
        ncolors=[min(subblockwidth,maximum(cols[i:blockwidth:nblock])) for i in 1:min(blockwidth,nblock)]
        endinds=cumsum(ncolors)
        startinds=[endinds[i]-ncolors[i]+1 for i in 1:min(blockwidth,nblock)]
        colors=[_cycle(startinds[blockcolors[i]]:endinds[blockcolors[i]],cols[i]) for i in 1:nblock]
        vcat(colors...)
    end

  end
end

end