text:
import MitosisStochasticDiffEq as MSDE
using Mitosis
using DiffEqNoiseProcess
using Test, Random
using LinearAlgebra
using Statistics
using StochasticDiffEq
using SparseArrays
using StaticArrays

# Test outer function
@testset "outer function tests" begin
  exA = rand(10,10)
  @test minimum(MSDE.outer_(exA) .!= 0)
  exB = Diagonal(exA)
  @test sum(MSDE.outer_(exB) .!= 0) == 10
end

"""
    forwardguiding(M, s, x, ps, Z) -> xs, ll
Forward sample a guided trajectory `xs` starting in `x` and compute it's
log-likelihood `ll` with innovations `Z = randn(length(s))`.
"""
function forwardguiding(plin, pest, s, (x, ll), ps, Z=randn(length(s)), noisetype=:scalar)
    # linear approximation of b and constant approximation of σ
    # with parameters B, β, and σ̃
    flinear(u,p,t) = p[1]*u .+ p[2]
    σlinear(u,p,t) = p[3]

    function llstep(x, r, t, P, noisetype)
      tmp = MSDE.outer_(g(x,pest,t)) - MSDE.outer_(σlinear(x,plin,t))
      dll = dot(f(x,pest,t) - flinear(x,plin,t), r) -0.5*tr(tmp*(inv(P) - MSDE.outer_(r)))
    end

    xs = typeof(x)[]
    d = length(x)
    for i in eachindex(s)[1:end-1]
        dt = s[i+1] - s[i]
        t = s[i]
        push!(xs, x)
        # ν = @view ps[:,i][1:d]
        # P = reshape(@view(ps[:,i][d+1:d+d*d]), d, d)
        # r = inv(P)*(ν .- x)
        ν = ps[i][1]
        P = ps[i][2]
        r = inv(P)*(ν .- x)

        ll += llstep(x, r, t, P, noisetype)*dt # accumulate log-likelihood

        if noisetype == :scalar
            noise = g(x,pest,t)*Z[i][] #sqrt(dt)*Z[i]
        elseif noisetype ==:diag
            noise = g(x,pest,t).*Z[i]
        elseif noisetype ==:nondiag
            noise = g(x,pest,t)*Z[:,i]
        else
            error("noisetype not understood.")
        end
        if x isa Number
           tmp = (MSDE.outer_(g(x,pest,t))*r*dt)[1]
       else
           tmp = MSDE.outer_(g(x,pest,t))*r*dt
       end
        x = x + f(x,pest,t)*dt + tmp + noise # evolution guided by observations

    end
    push!(xs, x)
    xs, ll
end

# define SDE function
f(u,p,t) = @. p[1]*u + p[2] - 1.5*sin(u*2pi)
g(u,p,t) = p[3] .- 0.2*(1 .-sin.(u))

@testset "IIP Guiding tests" begin
  # set true model parameters
  p = [-0.1,0.2,0.9]

  # set of linear parameters Eq.~(2.2)
  B, β, σ̃ = -0.1, 0.2, 1.3
  plin = [B, β, σ̃]
  pest = [-0.4, 0.5, 1.4] # initial guess of parameter to be estimated

  # time span
  tstart = 0.0
  tend = 1.0
  dt = 0.001
  trange = tstart:dt:tend

  # intial condition
  u0 = 1.1

  # forward kernel
  sdekernel = MSDE.SDEKernel(f,g,trange,pest)

  # initial values for ODE
  mynames = (:logscale, :μ, :Σ);
  myvalues = [0.0, 0.0, 10.0];
  NT = NamedTuple{mynames}(myvalues)

  # backward kernel
  kerneltilde = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, plin)
  message, backward = MSDE.backwardfilter(kerneltilde, NT)

  x0 = randn()
  ll0 = randn()

  (uend1, ll), (ts1, u1, noise1) = MSDE.forwardguiding(sdekernel, message, (x0, ll0); save_noise=true)

  dWs = (noise1[2:end]-noise1[1:end-1])
  ps = message.soldis
  solfw2, ll2 = forwardguiding(plin, pest, message.ts, (x0, ll0),ps,dWs)

  @test isapprox(getindex.(u1,1), solfw2, rtol=1e-12)
  @test isapprox(ll, ll2, rtol=1e-12)

  # multivariate tests with scalar random process
  dim = 7
  Random.seed!(1234)
  logscale = randn()
  μ = randn(dim)
  Σ = randn(dim,dim)
  myvalues = [logscale, μ, Σ];
  NT = NamedTuple{mynames}(myvalues)

  m = 1
  plin = [randn(dim,dim), randn(dim), randn(dim,m)] # B, β, σtil

  # set scalar random process
  t = tstart:dt:tend
  W = sqrt(dt)*randn(length(t))
  W1 = cumsum([zero(dt); W[1:end-1]])
  NG = NoiseGrid(t,W1)

  # backward kernel
  kerneltilde = MSDE.SDEKernel(Mitosis.AffineMap(plin[1], plin[2]), Mitosis.ConstantMap(plin[3]), trange, plin)
  message, backward = MSDE.backwardfilter(kerneltilde, NT)

  x0 = randn(dim)
  ll0 = randn()
  (uend1, ll), (ts1, u1, noise1) = MSDE.forwardguiding(sdekernel, message, (x0, ll0), EM(false), NG)

  ps = message.soldis
  solfw2, ll2 = forwardguiding(plin, pest, message.ts, (x0, ll0), ps, W)

  @test isapprox(u1, solfw2, rtol=1e-12)
  @test isapprox(ll, ll2, rtol=1e-12)

  # multivariate tests with diagonal noise random process
  dim = 2
  Random.seed!(12345)
  logscale = randn()
  μ = randn(dim)
  Σ = randn(dim,dim)
  myvalues = [logscale, μ, Σ];
  NT = NamedTuple{mynames}(myvalues)

  m = 2
  plin = [randn(dim,dim), randn(dim), randn(dim,m)] # B, β, σtil

  # backward kernel
  kerneltilde = MSDE.SDEKernel(Mitosis.AffineMap(plin[1], plin[2]), Mitosis.ConstantMap(plin[3]), trange, plin)
  message, backward = MSDE.backwardfilter(kerneltilde, NT)

  x0 = randn(dim)
  ll0 = randn()

  (uend1, ll), (ts1, u1, noise1) = MSDE.forwardguiding(sdekernel, message, (x0, ll0); save_noise=true)

  dWs = noise1[2:end]-noise1[1:end-1]

  ps = message.soldis
  solfw2, ll2 = forwardguiding(plin, pest, message.ts, (x0, ll0),ps,dWs,:diag)

  @test isapprox(u1, solfw2, rtol=1e-12)
  @test isapprox(ll, ll2, rtol=1e-12)
end

@testset "OOP Guiding tests" begin

  # set true model parameters
  p = [-0.1,0.2,0.9]

  # set of linear parameters Eq.~(2.2)
  B, β, σ̃ = -0.1, 0.2, 1.3
  plin = [B, β, σ̃]
  pest = [-0.4, 0.5, 1.4] # initial guess of parameter to be estimated

  # time span
  tstart = 0.0
  tend = 1.0
  dt = 0.001
  trange = tstart:dt:tend

  # intial condition
  u0 = 1.1

  # forward kernel
  sdekernel = MSDE.SDEKernel(f,g,trange,pest)

  # initial values for ODE
  mynames = (:logscale, :μ, :Σ);
  myvalues = [0.0, 0.0, 10.0];
  NT = NamedTuple{mynames}(myvalues)

  # backward kernel
  kerneltilde = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, plin)
  message, backward = MSDE.backwardfilter(kerneltilde, NT)

  x0 = randn()
  ll0 = randn()

  (uend1, ll), (ts1, u1, noise1) = MSDE.forwardguiding(sdekernel, message, (x0, ll0); save_noise=true, inplace=false)

  dWs = (noise1[2:end]-noise1[1:end-1])
  ps = message.soldis
  solfw2, ll2 = forwardguiding(plin, pest, message.ts, (x0, ll0),ps,dWs)

  @test isapprox(getindex.(u1,1), solfw2, rtol=1e-12)
  @test isapprox(ll, ll2, rtol=1e-12)

  # multivariate tests with scalar random process
  dim = 7
  Random.seed!(1234)
  logscale = randn()
  μ = randn(dim)
  Σ = randn(dim,dim)
  myvalues = [logscale, μ, Σ];
  NT = NamedTuple{mynames}(myvalues)

  m = 1
  plin = [randn(dim,dim), randn(dim), randn(dim,m)] # B, β, σtil

  # set scalar random process
  t = tstart:dt:tend
  W = sqrt(dt)*randn(length(t))
  W1 = cumsum([zero(dt); W[1:end-1]])
  NG = NoiseGrid(t,W1)

  # backward kernel
  kerneltilde = MSDE.SDEKernel(Mitosis.AffineMap(plin[1], plin[2]), Mitosis.ConstantMap(plin[3]), trange, plin)
  message, backward = MSDE.backwardfilter(kerneltilde, NT)

  x0 = randn(dim)
  ll0 = randn()
  (uend1, ll), (ts1, u1, noise1)  = MSDE.forwardguiding(sdekernel, message, (x0, ll0), EM(false), NG, inplace=false)

  ps = message.soldis
  solfw2, ll2 = forwardguiding(plin, pest, message.ts, (x0, ll0), ps, W)

  @test isapprox(u1, solfw2, rtol=1e-12)
  @test isapprox(ll, ll2, rtol=1e-12)

  # multivariate tests with diagonal noise random process
  dim = 2
  Random.seed!(12345)
  logscale = randn()
  μ = randn(dim)
  Σ = randn(dim,dim)
  myvalues = [logscale, μ, Σ];
  NT = NamedTuple{mynames}(myvalues)

  m = 2
  plin = [randn(dim,dim), randn(dim), randn(dim,m)] # B, β, σtil

  # backward kernel
  kerneltilde = MSDE.SDEKernel(Mitosis.AffineMap(plin[1], plin[2]), Mitosis.ConstantMap(plin[3]), trange, plin)
  message, backward = MSDE.backwardfilter(kerneltilde, NT)

  x0 = randn(dim)
  ll0 = randn()

  (uend1, ll), (ts1, u1, noise1)  = MSDE.forwardguiding(sdekernel, message, (x0, ll0); save_noise=true, inplace=false)

  Ws = Array(noise1)
  dWs = Ws[2:end]-Ws[1:end-1]

  ps = message.soldis
  solfw2, ll2 = forwardguiding(plin, pest, message.ts, (x0, ll0),ps,dWs,:diag)

  @test isapprox(u1, solfw2, rtol=1e-12)
  @test isapprox(ll, ll2, rtol=1e-12)
end

@testset "Adaptive Guiding tests" begin
  Random.seed!(12345)
  using StochasticDiffEq, DiffEqNoiseProcess
  # set true model parameters
  p = [-0.1,0.2,0.9]

  # set of linear parameters Eq.~(2.2)
  B, β, σ̃ = -0.1, 0.2, 1.3
  plin = [B, β, σ̃]
  pest = [-0.4, 0.5, 1.4] # initial guess of parameter to be estimated

  # time span
  tstart = 0.0
  tend = 1.0
  dt = 0.001
  trange = tstart:dt:tend

  # intial condition
  u0 = 1.1

  # forward kernel
  sdekernel = MSDE.SDEKernel(f,g,trange,pest)

  # initial values for ODE
  mynames = (:logscale, :μ, :Σ);
  myvalues = [0.0, 0.0, 10.0];
  NT = NamedTuple{mynames}(myvalues)

  # backward kernel
  kerneltilde = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, plin)
  message, backward = MSDE.backwardfilter(kerneltilde, NT)

  # define NoiseGrid
  brownian_values = cumsum([[zeros(2)];[sqrt(dt)*randn(2) for i in 1:length(trange)-1]])
  brownian_values2 = cumsum([[zeros(2)];[sqrt(dt)*randn(2) for i in 1:length(trange)-1]])
  W = NoiseGrid(collect(trange),brownian_values,brownian_values2)

  x0 = randn()
  ll0 = randn()

  (uend1, ll), (ts1, u1, noise1) = MSDE.forwardguiding(sdekernel, message, (x0, ll0), LambaEM(),
            W; dt=dt, isadaptive=false)
  (uend2, ll2), (ts2, u2, noise2) = MSDE.forwardguiding(sdekernel, message, (x0, ll0), LambaEM(),
            W; dt=dt, isadaptive=true)
  (uend3, ll3), (ts3, u3, noise3) = MSDE.forwardguiding(sdekernel, message, (x0, ll0), SOSRI(),
            W; dt=dt, isadaptive=true)

  @test isapprox(ll, ll2, rtol=1e-1)
  @test isapprox(ll, ll3, rtol=1e-1)
  @test isapprox(ll2, ll2, rtol=1e-1)

  @test length(ts1) == length(trange)
  @test length(ts2) < length(ts1)
  @test length(ts3) < length(ts1)

  @show length(ts1), length(ts2), length(ts3)
end

@testset "timechange Guiding tests" begin
  Random.seed!(12345)
  # set true model parameters
  p = [-0.1,0.2,0.9]

  # set of linear parameters Eq.~(2.2)
  B, β, σ̃ = -0.1, 0.2, 1.3
  plin = [B, β, σ̃]
  pest = [-0.4, 0.5, 1.4] # initial guess of parameter to be estimated

  # time span
  tstart = 0.0
  tend = 1.0
  dt = 0.001
  trange = tstart:dt:tend

  # intial condition
  u0 = 1.1

  # forward kernel
  sdekernel = MSDE.SDEKernel(f,g,trange,pest)

  # initial values for ODE
  mynames = (:logscale, :μ, :Σ);
  myvalues = [0.0, 0.0, 10.0];
  NT = NamedTuple{mynames}(myvalues)

  # backward kernel
  kerneltilde = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, plin)
  message, backward = MSDE.backwardfilter(kerneltilde, NT, apply_timechange=true)

  x0 = randn()
  ll0 = randn()

  (uend, ll), (ts, u, noise) = MSDE.forwardguiding(sdekernel, message, (x0, ll0);
    isadaptive=false)

  @test isapprox(ts, message.ts, rtol=1e-10)
  @test isapprox(ts, MSDE.timechange(trange), rtol=1e-10)
  @test length(ts) == length(trange)
end

@testset "Reuse of noise values tests" begin
  Random.seed!(12345)
  using StochasticDiffEq, DiffEqNoiseProcess
  # set true model parameters
  p = [-0.1,0.2,0.9]

  # set of linear parameters Eq.~(2.2)
  B, β, σ̃ = -0.1, 0.2, 1.3
  plin = [B, β, σ̃]
  pest = [-0.4, 0.5, 1.4] # initial guess of parameter to be estimated

  # time span
  tstart = 0.0
  tend = 1.0
  dt = 0.001
  trange = tstart:dt:tend

  # intial condition
  u0 = 1.1

  # forward kernel
  sdekernel = MSDE.SDEKernel(f,g,trange,pest)

  # initial values for ODE
  mynames = (:logscale, :μ, :Σ);
  myvalues = [0.0, 0.0, 10.0];
  NT = NamedTuple{mynames}(myvalues)

  # backward kernel
  kerneltilde = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, plin)
  message, backward = MSDE.backwardfilter(kerneltilde, NT)

  # define NoiseGrid
  brownian_values = cumsum([[zeros(2)];[sqrt(dt)*randn(2) for i in 1:length(trange)-1]])
  W = NoiseGrid(collect(trange),brownian_values)

  x0 = randn()
  ll0 = randn()

  # test two subsequent evaluations with same Brownian motion given by NoiseGrid
  (uend1, ll1), (ts1, u1, noise1) = MSDE.forwardguiding(sdekernel, message, (x0, ll0), EM(), W; dt=dt)
  (uend2, ll2), (ts2, u2, noise2) = MSDE.forwardguiding(sdekernel, message, (x0, ll0), EM(), W; dt=dt)

  @test isapprox(ll1, ll2, rtol=1e-14)
  @test isapprox(u1, u2, rtol=1e-14)
  @test isapprox(noise1, noise2, rtol=1e-14)
  @test isapprox(noise1, W.W, rtol=1e-14)

  # test pCN with \rho = 1
  (uend1, ll1), (ts1, u1, noise1) = MSDE.forwardguiding(sdekernel, message, (x0, ll0), EM(); dt=dt, save_noise=true)
  Z = pCN(ts1, noise1, 1.0)

  (uend2, ll2), (ts2, u2, noise2) = MSDE.forwardguiding(sdekernel, message, (x0, ll0), EM(), Z; dt=dt)

  @test isapprox(ll1, ll2, rtol=1e-14)
  @test isapprox(u1, u2, rtol=1e-14)
  @test isapprox(noise1, noise1, rtol=1e-14)
  @test isapprox(noise1, Z.W, rtol=1e-14)

  # test pCN with ρ = 0.2 (decrease dt for test)
  ρ = 0.2
  dt = 0.0001
  trange = tstart:dt:tend
  # backward kernel
  kerneltilde = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, plin)
  message, backward = MSDE.backwardfilter(kerneltilde, NT)
  (uend1, ll1), (ts1, u1, noise1) = MSDE.forwardguiding(sdekernel, message, (x0, ll0), EM(); dt=dt, save_noise=true)

  Z = pCN(ts1, noise1, ρ)

  (uend2, ll2), (ts2, u2, noise2) = MSDE.forwardguiding(sdekernel, message, (x0, ll0), EM(), Z; dt=dt)

  computedW(W, indx, dt) = (W[indx+1][1]-W[indx][1])/sqrt(dt) # likelihood part can be ignored
  dWnew = []
  dWold = []
  for i in 1:length(ts1[1:end-1])
    push!(dWnew,computedW(noise1,i,dt))
    push!(dWold,computedW(noise2,i,dt))
  end
  @show cor(dWnew,dWold)
  @test ≈(cor(dWnew,dWold),ρ,rtol=1e-1)
end

@testset "flag-constant and matrix-valued diffusivity tests" begin

  Random.seed!(12345)

  ## define model (two traits, as in phylo)
  d = 2
  bθ = [.5, 0.9]
  σ0 = [1.25, 1.35]
  θ0 = (bθ, σ0)

  θ = ([-1.0, -0.2], [0.1, 0.1])
  M = -0.05I + [-1.0 1.0; 1.0 -1.0]

  B(θ) = Diagonal(θ[1]) * M
  Σ(θ) = Diagonal(θ[2])

  f(u,θ,t) = tanh.(Diagonal(θ[1]) * M * u)
  g(u,θ,t) = θ[2]

  u0 = zeros(2)

  θlin = (B(θ), zeros(d), Σ(θ))

  # time span
  tstart = 0.0
  tend = 0.1
  dt = 0.001
  trange = tstart:dt:tend

  # forward kernels
  κ1 = MSDE.SDEKernel(f, g, trange, θ0)
  κ2 = MSDE.SDEKernel(f, g, trange, θ0, nothing, true)
  # backward kernel
  κ̃ = MSDE.SDEKernel(Mitosis.AffineMap(θlin[1], θlin[2]), Mitosis.ConstantMap(θlin[3]), trange, θlin)

  # forward sample
  xT, (ts, x, noise) = MSDE.sample(κ1, u0; save_noise=true)
  Z = pCN(ts, noise, 1)
  xT2, (_, x2, _) = MSDE.sample(κ2, u0, EM(false), Z)

  @test x ≈ x2
  @test xT ≈ xT2

  # backward filter
  logscale = randn()
  ν = randn(d)
  P = randn(d,d)
  gaussian = WGaussian{(:μ,:Σ,:c)}(ν, P, logscale)

  message, backward = MSDE.backwardfilter(κ̃, gaussian)

  # forward guiding
  κg1 = MSDE.SDEKernel(f, g, trange, θ)
  κg2 = MSDE.SDEKernel(f, g, trange, θ, nothing, true)

  x0 = randn(d)
  ll0 = randn()

  (uend1, ll1), (ts1, u1, noise1) = MSDE.forwardguiding(κg1, message, (x0, ll0); save_noise=true)
  Z = pCN(ts1, noise1, 1.0)
  (uend2, ll2), (ts2, u2, noise2) = MSDE.forwardguiding(κg2, message, (x0, ll0), EM(false), Z)

  @test ll1 ≈ ll2
  @test isapprox(u1, u2, rtol=1e-14)
  @test isapprox(noise1, noise2, rtol=1e-14)

  # check guiding with matrix-valued diffusion
  gmat(u,θ,t) = Diagonal(θ[2])
  kg3 = MSDE.SDEKernel(f, gmat, trange, θ, Σ(θ), true)

  # inplace=true
  Z = pCN(ts1, noise1, 1.0)
  (uend3, ll3), (ts3, u3, noise3) = MSDE.forwardguiding(kg3, message, (x0, ll0), EM(false), Z)
  @test ll1 ≈ ll3
  @test isapprox(u1, u3, rtol=1e-14)

  # inplace=false
  Z = pCN(ts1, noise1, 1.0)
  (uend3, ll3), (ts3, u3, noise3) = MSDE.forwardguiding(kg3, message, (x0, ll0), EM(false), Z, inplace=false)
  @test ll1 ≈ ll3
  @test isapprox(u1, u3, rtol=1e-14)
end

@testset "multivariate forward guiding tests" begin
  seed = 12345
  Random.seed!(seed)
  d = 2
  u0 = randn(2)
  θlin = (randn(d,d), randn(d), Diagonal([0.1, 0.1]))

  Σ(θ) = Diagonal(θ[2]) # just to generate the noise_rate_prototype

  f(u,p,t) = p[1]*u + p[2]
  gvec(u,p,t) = diag(p[3])
  g(u,p,t) = p[3]
  function gstepvec!(dx, _, u, p, t, dw, _)
    dx .+= diag(p[3]).*dw
  end

  function gstep!(dx, _, u, p, t, dw, _)
    dx .+= p[3]*dw
  end

  # Define a sparse matrix by making a dense matrix and setting some values as not zero
  A = zeros(2,2)
  A[1,1] = 1
  A[2,2] = 1
  A = sparse(A)

  # time range
  tstart = 0.0
  tend = 1.0
  dt = 0.02
  trange = tstart:dt:tend

  # define kernels
  k1 = MSDE.SDEKernel(f,gvec,trange,θlin)
  k2 = MSDE.SDEKernel(f,g,trange,θlin,Σ(θlin))
  k3 = MSDE.SDEKernel(Mitosis.AffineMap(θlin[1], θlin[2]), Mitosis.ConstantMap(θlin[3]), trange, θlin, Σ(θlin))

  solend, (_, sol, _) = MSDE.sample(k1, u0, EM(false), save_noise=true)
  v = solend
  c = randn()
  Pmat = randn(d,d)

  message1, backward1 = MSDE.backwardfilter(k1, WGaussian{(:μ, :Σ, :c)}(v, Pmat, c))
  message2, backward2 = MSDE.backwardfilter(k2, WGaussian{(:μ, :Σ, :c)}(v, Pmat, c))
  message3, backward3 = MSDE.backwardfilter(k2, WGaussian{(:μ, :Σ, :c)}(v, Pmat, c))
  message4, backward4 = MSDE.backwardfilter(k3, WGaussian{(:μ, :Σ, :c)}(SVector{length(v)}(v), SMatrix{2,2,eltype(Pmat)}(Pmat), c))

  @test message1.soldis == message2.soldis
  @test message1.soldis == message3.soldis
  @test getindex.(message1.soldis,1) ≈ getindex.(message4.soldis,1) rtol=1e-14
  @test getindex.(message1.soldis,2) ≈ getindex.(message4.soldis,2) rtol=1e-14
  @test getindex.(message1.soldis,3) ≈ getindex.(message4.soldis,3) rtol=1e-14

  ll0 = randn()

  @testset "StochasticDiffEq EM() solver" begin
    (uend1, ll1), (ts1, u1, noise1) = MSDE.forwardguiding(k1, message1, (u0, ll0), EM(false); save_noise=true)
    Z = pCN(ts1, noise1, 1.0)
    (uend2, ll2), (ts2, u2, noise2) = MSDE.forwardguiding(k2, message2, (u0, ll0), EM(false), Z; save_noise=true)
    Z = pCN(ts1, noise1, 1.0)
    (uend3, ll3), (ts3, u3, noise3) = MSDE.forwardguiding(k3, message3, (u0, ll0), EM(false), Z; save_noise=true)
    Z = pCN(ts1, noise1, 1.0)
    (uend4, ll4), (ts4, u4, noise4) = MSDE.forwardguiding(k3, message4, (u0, ll0), EM(false), Z; save_noise=true)

    @test u1 ≈ u2 rtol=1e-14
    @test u1 ≈ u3 rtol=1e-14
    @test u1 ≈ u4 rtol=1e-14
    @test ll1 == ll2
    @test ll1 == ll3
    @test ll1 == ll4
  end

  @testset "internal solver" begin
    @testset "without passing a noise" begin
      Random.seed!(seed)
      (uend1, ll1), (ts1, u1, noise1) = MSDE.forwardguiding(k1, message1, (u0, ll0), MSDE.EulerMaruyama!(), inplace=false)
      Random.seed!(seed)
      (uend2, ll2), (ts2, u2, noise2) = MSDE.forwardguiding(k2, message2, (u0, ll0), MSDE.EulerMaruyama!(), inplace=false)
      Random.seed!(seed)
      (uend3, ll3), (ts3, u3, noise3) = MSDE.forwardguiding(k3, message3, (u0, ll0), MSDE.EulerMaruyama!(), inplace=false)
      Random.seed!(seed)
      (uend4, ll4), (ts4, u4, noise4) = MSDE.forwardguiding(k3, message4, (u0, ll0), MSDE.EulerMaruyama!(), inplace=false)
      @test isapprox(uend1,uend2,rtol=1e-14)
      @test isapprox(uend1,uend3,rtol=1e-14)
      @test isapprox(uend1,uend4,rtol=1e-14)
      @test ll1 == ll2
      @test ll1 == ll3
      @test ll1 == ll4
    end

    @testset "passing a noise grid" begin
      Ws = cumsum([[zero(u0)];[sqrt(trange[i+1]-ti)*randn(size(u0))
                for (i,ti) in enumerate(trange[1:end-1])]])
      NG = NoiseGrid(trange,Ws)

      (uendEM, llEM), (tsEM, uEM, noiseEM) = MSDE.forwardguiding(k3, message4, (u0, ll0), EM(false), NG, inplace=false)
      (uend1, ll1), (ts1, u1, noise1) = MSDE.forwardguiding(k1, message1, (u0, ll0), MSDE.EulerMaruyama!(), NG, inplace=false)
      (uend2, ll2), (ts2, u2, noise2) = MSDE.forwardguiding(k2, message2, (u0, ll0), MSDE.EulerMaruyama!(), NG, inplace=false)
      (uend3, ll3), (ts3, u3, noise3) = MSDE.forwardguiding(k3, message3, (u0, ll0), MSDE.EulerMaruyama!(), NG, inplace=false)
      (uend4, ll4), (ts4, u4, noise4) = MSDE.forwardguiding(k3, message4, (u0, ll0), MSDE.EulerMaruyama!(), NG, inplace=false)

      @test u1 ≈ uEM rtol=1e-12
      @test llEM ≈ ll1 rtol=1e-12

      @test isapprox(uend1,uend2,rtol=1e-14)
      @test isapprox(uend1,uend3,rtol=1e-14)
      @test isapprox(uend1,uend4,rtol=1e-14)
      @test ll1 == ll2
      @test ll1 == ll3
      @test ll1 == ll4
    end

    @testset "passing the noise values" begin
      # pass noise process and compare with EM()
      Ws = cumsum([[zero(u0)];[sqrt(trange[i+1]-ti)*randn(size(u0))
              for (i,ti) in enumerate(trange[1:end-1])]])
      NG = NoiseGrid(trange,Ws)

      (uendEM, llEM), (tsEM, uEM, noiseEM) = MSDE.forwardguiding(k3, message4, (u0, ll0), EM(false), NG, inplace=false)
      (uend1, ll1), (ts1, u1, noise1)  = MSDE.forwardguiding(k1, message1, (u0, ll0), MSDE.EulerMaruyama!(), Ws, inplace=false)
      (uend2, ll2), (ts2, u2, noise2) = MSDE.forwardguiding(k2, message2, (u0, ll0), MSDE.EulerMaruyama!(), Ws, inplace=false)
      (uend3, ll3), (ts3, u3, noise3) = MSDE.forwardguiding(k3, message3, (u0, ll0), MSDE.EulerMaruyama!(), Ws, inplace=false)
      (uend4, ll4), (ts4, u4, noise4) = MSDE.forwardguiding(k3, message4, (u0, ll0), MSDE.EulerMaruyama!(), Ws, inplace=false)

      @test u1 ≈ uEM rtol=1e-12
      @test llEM ≈ ll1 rtol=1e-12

      @test isapprox(uend1,uend2,rtol=1e-14)
      @test isapprox(uend1,uend3,rtol=1e-14)
      @test isapprox(uend1,uend4,rtol=1e-14)
      @test ll1 == ll2
      @test ll1 == ll3
      @test ll1 == ll4
    end
  end
end


text:
# https://fenicsproject.org/docs/dolfin/latest/python/demos/cahn-hilliard/demo_cahn-hilliard.py.html

module cahn_hilliard

using FenicsPy

py"""
import random
from dolfin import UserExpression, MPI, NonlinearProblem, assemble

# Class representing the intial conditions
class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        random.seed(2 + MPI.rank(MPI.comm_world))
        super().__init__(**kwargs)
    def eval(self, values, x):
        values[0] = 0.63 + 0.02*(0.5 - random.random())
        values[1] = 0.0
    def value_shape(self):
        return (2,)

# Class for interfacing with the Newton solver
class CahnHilliardEquation(NonlinearProblem):
    def __init__(self, a, L):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=A)
"""

# Model parameters
lmbda  = 1.0e-02  # surface parameter
dt     = 5.0e-06  # time step
theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson

# Form compiler options
parameters["form_compiler"]["optimize"]     = true
parameters["form_compiler"]["cpp_optimize"] = true

# Create mesh and build function space
mesh = UnitSquareMesh.create(96, 96, CellType.Type.quadrilateral)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)

# Define trial and test functions
du    = TrialFunction(ME)
q, v  = TestFunctions(ME)

# Define functions
u   = FeFunction(ME)  # current solution
u0  = FeFunction(ME)  # solution from previous converged step

# Split mixed functions
c,  mu  = split(u)
c0, mu0 = split(u0)
dc, dmu = split(du)

# Create intial conditions and interpolate
u_init = Expression(py"InitialConditions"(degree=1))
u.interpolate(u_init)
u0.interpolate(u_init)

# Compute the chemical potential df/dc
c = variable(c)
f    = 100*c^2*(1-c)^2
dfdc = diff(f, c)

# mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu

# Weak statement of the equations
L0 = c*q*dx - c0*q*dx + dt*dot(grad(mu_mid), grad(q))*dx
L1 = mu*v*dx - dfdc*v*dx - lmbda*dot(grad(c), grad(v))*dx
L = L0 + L1

# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)

# Create nonlinear problem and Newton solver
problem = NonlinearProblem(py"CahnHilliardEquation"(a.pyobject, L.pyobject))
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

# Output file
file = File("cahn_hilliard/output.pvd", "compressed")

# Step in time
t = 0.0
T = 50*dt
while t < T
    global t += dt
    solver.solve(problem, u.vector())
    u0.assign(u)
    file << (u.split()[1], t)
end

end # module cahn_hilliard

ext: jl


text:
################################################################################
### CUSTOM DATA STRUCTURES
#=
Custom data structures are defined here
(structures are immutable by default; if modification is required, use "mutable struct")
=#

################################################################################
### PRIMARY INDEX SETS
struct ArcDataStruct
    ID::Array
    node_s::Dict
    node_r::Dict
    cap::Dict
    len::Dict
end

struct NodeDataStruct
    ID::Array
    alias::Dict
    lon::Dict
    lat::Dict
end

struct ProductDataStruct
    ID::Array
    alias::Dict
    transport_cost::Dict
end

struct DemandDataStruct
    ID::Array
    node::Dict
    prod::Dict
    bid::Dict
    cap::Dict
end

struct SupplyDataStruct
    ID::Array
    node::Dict
    prod::Dict
    bid::Dict
    cap::Dict
end

struct TechDataStruct
    ID::Array
    Outputs::Dict
    Inputs::Dict
    OutputYields::Dict
    InputYields::Dict
    InputRef::Dict
    bid::Dict
    cap::Dict
    alias::Dict
end

struct TechsiteDataStruct
    ID::Array
    node::Dict
    tech::Dict
end

################################################################################
### SECONDARY SETS (SUBSETS, OR INTERSECTIONS OF TWO OR MORE PRIMARY INDICES)
struct SetStruct
    Ain::Dict
    Aout::Dict
    TPQ::Array
    NQT::Dict
    NPt::Dict
    NQt::Dict
end

################################################################################
### PARAMETERS
struct ParStruct
    dMAX::Dict
    sMAX::Dict
    α::Dict
    ξgenMAX::Dict
    ξconMAX::Dict
    fMAX::Dict
end

################################################################################
### MODEL OUTPUTS
struct OutputStruct
    z::Float64
    si::JuMP.Containers.DenseAxisArray
    dj::JuMP.Containers.DenseAxisArray
    snp::JuMP.Containers.DenseAxisArray
    dnp::JuMP.Containers.DenseAxisArray
    f::JuMP.Containers.DenseAxisArray
    ξcon::JuMP.Containers.DenseAxisArray
    ξgen::JuMP.Containers.DenseAxisArray
    πNP::JuMP.Containers.DenseAxisArray
    πA::Dict
    πT::Dict
    Φd::Dict
    Φs::Dict
    Φf::Dict
    Φξ::Dict
end

struct StatStruct
    NumVars::Float64
    TotalIneqCons::Float64
    TotalEqCons::Float64
    NumVarBounds::Float64
    ModelIneqCons::Float64
    ModelEqCons::Float64
end


text:
include("find.jl")
include("fit.jl")

max_stars_repo_name: HolyLab/ExtractPSF

max_forks_repo_name: HolyLab/ExtractPSF


text:
print_all(Lexers.FortranLexer, "fortran")

tokentest(
    Lexers.FortranLexer,
    "! a single line comment.",
    COMMENT => "! a single line comment.",
)

tokentest(
    Lexers.FortranLexer,
    "x = 1 ! a single line comment.",
    NAME => "x",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    NUMBER_INTEGER => "1",
    TEXT => " ",
    COMMENT => "! a single line comment.",
)

tokentest(
    Lexers.FortranLexer,
    "program test\n    print *, 'hello world'",
    KEYWORD => "program",
    TEXT => " ",
    NAME => "test",
    TEXT => "\n    ",
    KEYWORD => "print",
    TEXT => " ",
    OPERATOR => "*",
    PUNCTUATION => ",",
    TEXT => " ",
    STRING_SINGLE => "'hello world'",
)

tokentest(
    Lexers.FortranLexer,
    "print \"(F6.3)\", 4.32",
    KEYWORD => "print",
    TEXT => " ",
    STRING_DOUBLE => "\"(F6.3)\"",
    PUNCTUATION => ",",
    TEXT => " ",
    NUMBER_FLOAT => "4.32",
)

tokentest(
    Lexers.FortranLexer,
    "# a preprocessor comment",
    COMMENT_PREPROC => "# a preprocessor comment",
)

tokentest(
    Lexers.FortranLexer,
    ".true. == .FALSE. .and. .false.",
    NAME_BUILTIN => ".true.",
    TEXT => " ",
    OPERATOR => "==",
    TEXT => " ",
    NAME_BUILTIN => ".FALSE.",
    TEXT => " ",
    OPERATOR_WORD => ".and.",
    TEXT => " ",
    NAME_BUILTIN => ".false.",
)

tokentest(
    Lexers.FortranLexer,
    "TYPE string80\n    INTEGER length",
    KEYWORD => "TYPE",
    TEXT => " ",
    NAME => "string80",
    TEXT => "\n    ",
    KEYWORD_TYPE => "INTEGER",
    TEXT => " ",
    NAME => "length",
)

tokentest(
    Lexers.FortranLexer,
    "w = 5/v + a(1:5, 5)",
    NAME => "w",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    NUMBER_INTEGER => "5",
    OPERATOR => "/",
    NAME => "v",
    TEXT => " ",
    OPERATOR => "+",
    TEXT => " ",
    NAME => "a",
    PUNCTUATION => "(",
    NUMBER_INTEGER => "1",
    PUNCTUATION => ":",
    NUMBER_INTEGER => "5",
    PUNCTUATION => ",",
    TEXT => " ",
    NUMBER_INTEGER => "5",
    PUNCTUATION => ")",
)


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Parallel_Computing")

# # Parallel Algorithms: Thinking in Parallel
#
# Now that we're starting to see the challenges of parallelism, it's worth taking
# a step back and examining how we might go about designing parallel algorithms.
#
# This is adapted from a [workshop paper](http://jiahao.github.io/parallel-prefix/) by Jiahao Chen and
# Alan Edelman entitled "Parallel Prefix Polymorphism Permits Parallelization, Presentation & Proof" and
# will appear in the proceedings of the [First Workshop for High Performance Technical Computing in Dynamic
# Languages](http://jiahao.github.io/hptcdl-sc14/), held in conjunction with [SC14: The International Conference on High Performance Computing, Networking, Storage and Analysis](http://sc14.supercomputing.org/)

#-

using Compose, Gadfly

# # `reduce()`
#
# Reduction applies a binary operator to a vector repeatedly to return a scalar. Thus + becomes sum, and * becomes prod.
#
# It is considered a basic parallel computing primitive.

reduce(+, 1:8), sum(1:8)  # triangular numbers

#-

reduce(*, 1:8), prod(1:8) # factorials

#-

boring(a,b)=a
@show reduce(boring, 1:8)
boring2(a,b)=b
@show reduce(boring2, 1:8)

# You can also use reduce to compute Fibonacci numbers using their recurrences.
#
# $$\begin{pmatrix} f_2 \\f_1 \end{pmatrix} = \begin{pmatrix} f_1 + f_0 \\ f_1 \end{pmatrix}
# = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} f_1 \\ f_0 \end{pmatrix} $$
#
# $$\begin{pmatrix} f_{n+1} \\ f_n \end{pmatrix} = \dots
# = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n \begin{pmatrix} f_1 \\ f_0 \end{pmatrix} $$
#
# From this, you can show that
#
# $$\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n  = \begin{pmatrix} f_{n+1} & f_n \\ f_n & f_{n-1} \end{pmatrix} $$
#
# (this applies reduce to the same argument over and over again -- there are of course other ways)

M=[1 1; 1 0]
reduce(*,fill(M,3))
prod(fill(M,3))

#-

n= 50 # Try manipulating n to pick different values (try between 0-100)
prod(fill(big.(M),n))

#-

fib(j)=reduce(*, fill(M,j))
fib.([4,7])

# You can solve recurrences of any complexity using `reduce`. For example, `reduce` can compute a Hadamard matrix from its definition in terms of its submatrices:
#
# $$H_2 = \begin{pmatrix} H_1 & H_1 \\ H_1 & -H_1 \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} \otimes H_1$$
#
# and so on.
#
# (Note: this is just using reduce to compute a matrix power.
# One can think of alternative ways for sure.)

## If A is m x n
## If B is p x q
## then kron(A,B) is mp x nq and has all the elements of A times all of the elements of B

#-

A=[1 2;3 4]
B=[10 100; 1 -10]
⊗(A,B)=kron(A,B)

M=[ 1 1;1 -1]
H=⊗(⊗(⊗(M,M),M),M)

#-

H'H

#-

Hadamard(n)=reduce(⊗, fill(M,n))
H=Hadamard(4)

#-

using LinearAlgebra
cumsum(1:8)  # It is useful to know that cumsum is a linear operator
## You can use power method! Below is the underlying matrix
A=tril(ones(Int,8,8))

#-

[A*(1:8),cumsum(1:8)]

# # `prefix`

#-

# Having discussed `reduce`, we are now ready for the idea behind prefix sum.
# Prefix or scan is long considered an important parallel
# primitive as well.
#
# Suppose you wanted to compute the partial sums of a vector, i.e. given
# `y[1:n]`, we want to overwrite the vector `y` with the vector of partial sums
#
# ```julia
# new_y[1] = y[1]
# new_y[2] = y[1] + y[2]
# new_y[3] = y[1] + y[2] + y[3]
# ...
# ```
#
# At first blush, it seems impossible to parallelize this, since
#
# ```julia
# new_y[1] = y[1]
# new_y[2] = new_y[1] + y[2]
# new_y[3] = new_y[2] + y[3]
# ...
# ```
#
# which appears to be an intrinsically serial process. As written with a `+`
# operator, this is `cumsum` — but note that it can generalize to any operation.

function prefix_serial!(y, ⊕)
    for i=2:length(y)
        y[i] = y[i-1] ⊕ y[i]
    end
    y
end

#-

prefix_serial!([1:8;],+)

#-

cumsum(1:8)

#-

prefix_serial!([1:8;], *)

#-

cumprod(1:8)

# However, it turns out that because these operations are associative, we can regroup the _order_ of how these sums or products are carried out. (This of course extends to other associative operations, too.) Another ordering of 8 associative operations is provided by `prefix8!`:

## Magic :)
function prefix8!(y, ⊕)
    length(y)==8 || error("length 8 only")
    for i in (2,4,6,8); y[i] = y[i-1] ⊕ y[i]; end
    for i in (  4,  8); y[i] = y[i-2] ⊕ y[i]; end
    for i in (      8); y[i] = y[i-4] ⊕ y[i]; end
    for i in (    6  ); y[i] = y[i-2] ⊕ y[i]; end
    for i in ( 3,5,7 ); y[i] = y[i-1] ⊕ y[i]; end
    y
end

#-

prefix8!([1:8;], +) == cumsum(1:8)

# In fact, this can generalize beyond just length-8 arrays:

## More magic
function prefix!(y, ⊕)
    l=length(y)
    k=ceil(Int, log2(l))
    @inbounds for j=1:k, i=2^j:2^j:min(l, 2^k)              #"reduce"
        y[i] = y[i-2^(j-1)] ⊕ y[i]
    end
    @inbounds for j=(k-1):-1:1, i=3*2^(j-1):2^j:min(l, 2^k) #"expand"
        y[i] = y[i-2^(j-1)] ⊕ y[i]
    end
    y
end

# -

A = rand(0:9, 123)
prefix!(copy(A), *) == cumprod(A)

# ## What is this magic?

#-

# We can visualize the operations with a little bit of trickery. In Julia, arrays are simply types that expose the array protocol. In particular, they need to implement  methods for the generic functions `length`, `getindex` and `setindex!`. The last two are used in indexing operations, since statements
#
#     y = A[1]
#     A[3] = y
#
# get desugared to
#
#     y = getindex(A, 1)
#     setindex!(A, y, 3)
#
# respectively.
#
# We can trace through the iterable by introduce a dummy array type, `AccessArray`, which records every access to `getindex` and `setindex!`.
#
# Specifically:
#
# - `length(A::AccessArray)` returns the length of the array it wraps
# - `getindex(A::AccessArray, i)` records read access to the index `i` in the `A.read` field and then actually retuns the value in the array it wraps.
# - `setindex!(A::AccessArray, x, i)` records write access to the index `i`. The `A.history` field is appended with a new tuple consisting of the current `A.read` field and the index `i`, and then it performs the assignment.
#
# The way `AccessArray` works, it assumes an association between a single `setindex!` call and and all the preceding `getindex` calls since the previous `setindex!` call, which is sufficient for the purposes of tracing through prefix calls.

mutable struct AccessArray{T,N,A}
    data :: A
    read :: Vector{Int}
    history :: Vector{Tuple{Vector{Int},Int}}
end
AccessArray(A) = AccessArray{eltype(A), ndims(A), typeof(A)}(A, Int[], Int[])

Base.length(A::AccessArray) = length(A.data)

function Base.getindex(A::AccessArray, i::Int)
    push!(A.read, i)
    A.data[i]
end

function Base.setindex!(A::AccessArray, x, i::Int)
    push!(A.history, (A.read, i))
    A.read = Int[]
    A.data[i] = x
end

#-

M = AccessArray(rand(8))

#-

M[7] = M[3] + M[2]

#-

M.history

# So now we can trace the access pattern when calling `prefix8`!

A=prefix8!(AccessArray(rand(8)),+)

#-

A.history

# Now let's visualize this! Each entry in `A.history` is rendered by a gate object:

using Compose: circle, mm

#-

struct Gate{I,O}
    ins :: I
    outs :: O
end

import Gadfly.render

function render(G::Gate, x₁, y₁, y₀; rᵢ=0.1, rₒ=0.25)
    ipoints = [(i, y₀+rᵢ) for i in G.ins]
    opoints = [(i, y₀+0.5) for i in G.outs]
    igates  = [circle(i..., rᵢ) for i in ipoints]
    ogates  = [circle(i..., rₒ) for i in opoints]
    lines = [line([i, j]) for i in ipoints, j in opoints]
    compose(context(units=UnitBox(0.5,0,x₁,y₁+1)),
    compose(context(), stroke(colorant"black"), fill(colorant"white"),
            igates..., ogates...),
    compose(context(), linewidth(0.3mm), stroke(colorant"black"),
            lines...))
end

A=Gate([1,2],2)
render(A,2,0,0)

# Now we render the whole algorithm. We have to scan through the trace twice; the first time merely calculates the maximum depth that needs to be drawn and the second time actually generates the objects.

function render(A::AccessArray)
    #Scan to find maximum depth
    olast = depth = 0
    for y in A.history
        (any(y[1] .≤ olast)) && (depth += 1)
        olast = maximum(y[2])
    end
    maxdepth = depth

    olast = depth = 0
    C = []
    for y in A.history
        (any(y[1] .≤ olast)) && (depth += 1)
        push!(C, render(Gate(y...), length(A), maxdepth, depth))
        olast = maximum(y[2])
    end

    push!(C, compose(context(units=UnitBox(0.5,0,length(A),1)),
      [line([(i,0), (i,1)]) for i=1:length(A)]...,
    linewidth(0.1mm), stroke(colorant"grey")))
    compose(context(), C...)
end

#-

render(prefix!(AccessArray(zeros(8)), +))

# Now we can see that `prefix!` rearranges the operations to form two spanning trees:

render(prefix!(AccessArray(zeros(120)),+))

# as contrasted with the serial code:

render(prefix_serial!(AccessArray(zeros(8)),+))

#-

npp=90 # Try manipulating npp to choose different values between 1 and 180
render(prefix!(AccessArray(zeros(npp)),+))

#-

npp=90 # And manipulate it again, this time with the serial algorithm
render(prefix_serial!(AccessArray(zeros(npp)),+))

# # Now exploit the parallelism in the _algorithm_ to use a parallel _implementation_

using .Threads
function prefix_threads!(y, ⊕)
    l=length(y)
    k=ceil(Int, log2(l))
    for j=1:k
        @threads for i=2^j:2^j:min(l, 2^k)       #"reduce"
            @inbounds y[i] = y[i-2^(j-1)] ⊕ y[i]
        end
    end
    for j=(k-1):-1:1
        @threads for i=3*2^(j-1):2^j:min(l, 2^k) #"expand"
            @inbounds y[i] = y[i-2^(j-1)] ⊕ y[i]
        end
    end
    y
end

A = rand(500_000);

using BenchmarkTools
@btime prefix_serial!($(copy(A)), +);
@btime prefix!($(copy(A)), +);
@btime prefix_threads!($(copy(A)), +);

prefix_threads!(copy(A), +) == prefix!(copy(A), +) ≈ cumsum(A)

# # Thinking in parallel
#
# Notice how we didn't need to contort ourselves in making our algorithm
# work with `@threads`. We really did _just_ take a `@threads` on it and it
# just worked. It was both accurate _and_ fast.
#
# Coming up with rearrangements that make your particular algorithm parallel
# friendly isn't always easy, but when possible it makes everything else
# just fall out naturally.
#
# Finally, note that there can be clever ways to visualize algorithms as sanity checks.

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
using PartitionedArrays

function test_p_timers(parts)

  t = PTimer(parts)

  tic!(t;barrier=true)
  toc!(t,"Phase 1")
  toc!(t,"Phase 3")
  sleep(0.2)
  toc!(t,"Phase 143")
  tic!(t;barrier=true)
  sleep(0.4)
  toc!(t,"Matrix Assembly")

  display(t.data)

  map_main(t.data) do data
    open("times.txt","w") do io
      println(io,data)
    end
  end

  display(t)

  t = PTimer(parts,verbose=true)

  tic!(t;barrier=true)
  toc!(t,"Phase 1")
  toc!(t,"Phase 3")
  sleep(0.2)
  toc!(t,"Phase 143")
  tic!(t;barrier=true)
  sleep(0.4)
  toc!(t,"Matrix Assembly")

  display(t)

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

include("model_Alias.jl")
include("model_AliasPath.jl")
include("model_AliasPattern.jl")
include("model_BasicDependency.jl")
include("model_CloudError.jl")
include("model_DebugSetting.jl")
include("model_Dependency.jl")
include("model_Deployment.jl")
include("model_DeploymentExportResult.jl")
include("model_DeploymentExtended.jl")
include("model_DeploymentExtendedFilter.jl")
include("model_DeploymentListResult.jl")
include("model_DeploymentOperation.jl")
include("model_DeploymentOperationProperties.jl")
include("model_DeploymentOperationsListResult.jl")
include("model_DeploymentProperties.jl")
include("model_DeploymentPropertiesExtended.jl")
include("model_DeploymentValidateResult.jl")
include("model_DeploymentWhatIf.jl")
include("model_DeploymentWhatIfProperties.jl")
include("model_DeploymentWhatIfSettings.jl")
include("model_ErrorAdditionalInfo.jl")
include("model_ErrorResponse.jl")
include("model_ExportTemplateRequest.jl")
include("model_GenericResource.jl")
include("model_GenericResourceExpanded.jl")
include("model_GenericResourceFilter.jl")
include("model_HttpMessage.jl")
include("model_Identity.jl")
include("model_IdentityUserAssignedIdentities.jl")
include("model_OnErrorDeployment.jl")
include("model_OnErrorDeploymentExtended.jl")
include("model_Operation.jl")
include("model_OperationDisplay.jl")
include("model_OperationListResult.jl")
include("model_ParametersLink.jl")
include("model_Plan.jl")
include("model_Provider.jl")
include("model_ProviderListResult.jl")
include("model_ProviderResourceType.jl")
include("model_Resource.jl")
include("model_ResourceGroup.jl")
include("model_ResourceGroupExportResult.jl")
include("model_ResourceGroupFilter.jl")
include("model_ResourceGroupListResult.jl")
include("model_ResourceGroupPatchable.jl")
include("model_ResourceGroupProperties.jl")
include("model_ResourceListResult.jl")
include("model_ResourceProviderOperationDisplayProperties.jl")
include("model_ResourceReference.jl")
include("model_ResourcesMoveInfo.jl")
include("model_ScopedDeployment.jl")
include("model_ScopedDeploymentWhatIf.jl")
include("model_Sku.jl")
include("model_StatusMessage.jl")
include("model_SubResource.jl")
include("model_TagCount.jl")
include("model_TagDetails.jl")
include("model_TagValue.jl")
include("model_Tags.jl")
include("model_TagsListResult.jl")
include("model_TagsPatchResource.jl")
include("model_TagsResource.jl")
include("model_TargetResource.jl")
include("model_TemplateHashResult.jl")
include("model_TemplateLink.jl")
include("model_WhatIfChange.jl")
include("model_WhatIfOperationProperties.jl")
include("model_WhatIfOperationResult.jl")
include("model_WhatIfPropertyChange.jl")


text:
## The next programs shows how to use a if conditionals in advance way. Also shows how to invert them.

# Ags = it's the number of conditionals used in a program

# Build conditions such as:
# a == b;
# a == b && c == d               a == b || c == d;
# a == b && c == d && e == f     a == b || c == d || e == f ....;
# e.g.

# <QUANTITY> It's the number of varibales used in the computation. It must be major or equals to 2 and times 2.
# <return> Returns the all constructions with the values.
# <note> each variables has values true or false.

#Generate conditionals with ==
function generateEqualsConditionals(QUANTITY)
    if(!(QUANTITY%2==0))
        print("It's not possible compute\n")
    else
        print("Computing the conditionals\n")
    end
end

#Genreate conditionals with !=
function generateNOTEqualsConditionals(QUANTITY)
    if(!(QUANTITY%2==0))
        return
    else
        print("Todo el codigo")
    end
end

#Generate contditionals with == and !=
function generateMIXTUREConditionals(QUANTITY)
    if(!(QUANTITY%2==0))
        return
    else
        print("Todo el codigo")
    end
end

function printHeader()
    print("\n=============================================================================\tSTART\n")
end

function printFooter()
    print("\n========================== Programed by Madai Arteaga ========================\tEND\n")
end

##########################
## Executing main program

input_arg = ARGS[1]
printHeader()
print("This code is going to compute by $input_arg elements\n")
generateEqualsConditionals(parse(UInt8,input_arg))

printFooter()

max_stars_repo_name: MadaiA/JuliaProjects

max_issues_repo_name: MadaiA/JuliaProjects

max_forks_repo_name: MadaiA/JuliaProjects


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

# Compute X <- a X + b Y.
function matfun_axpby!(X,a,b,Y)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:m
        @inbounds for j=1:n
            if (b isa ValueOne)
                X[i,j]+=Y[i,j]
            else
                X[i,j]+=b*Y[i,j]
            end
        end
    end
end

@inline function exp_bbcs_m10(A)
    T=promote_type(eltype(A),ComplexF64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_bbcs_m10!(A_copy)
end

@inline function exp_bbcs_m10!(A)
    T=promote_type(eltype(A),ComplexF64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=7
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: Bb2 Ba2 B2 Bb3 B3 Ba5_4 B4 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 B8 B9 B10 B11 T2k13
    # Computing Bb2 = x*I+x*A
    coeff1=0.0 + 0.0im
    coeff2=0.03125 + 0.0im
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.0 + 0.0im
    coeff2=0.03125 + 0.0im
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Bb3 = x*A+x*B2
    coeff1=0.03125 + 0.0im
    coeff2=0.0 + 0.0im
    memslots2 .= coeff1.*memslots1 .+ coeff2.*memslots4
    # Computing B3 with operation: mult
    mul!(memslots3,memslots4,memslots2)
    # Deallocating Bb3 in slot 2
    # Computing Ba5_4 = x*A+x*B2+x*B3
    coeff1=0.00375 + 0.0im
    coeff2=0.0 - 0.008774760968797039im
    coeff3=-0.0009784845352378095 + 0.0im
    memslots2 .= coeff1.*memslots1 .+ coeff2.*memslots4 .+ coeff3.*memslots3
    # Computing B4 with operation: mult
    mul!(memslots5,memslots3,memslots3)
    # Computing Bb5 = x*B2+x*B3+x*B4
    coeff1=0.0 - 0.12395369585828313im
    coeff2=-0.011202694841085593 + 0.0im
    coeff3=-0.0 - 1.2367240538259895e-5im
    memslots6 .= coeff1.*memslots4 .+ coeff2.*memslots3 .+ coeff3.*memslots5
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots6)
    # Deallocating Ba5_4 in slot 2
    # Deallocating Bb5 in slot 6
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.3420232802536553 + 0.0im
    coeff2=0.0 - 0.008912493113538797im
    coeff3=0.047347067331271094 + 0.0im
    coeff4=0.0 - 0.022186600635366212im
    coeff5=-9.74758985615379e-6 + 0.0im
    coeff6=1.0 + 0.0im
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots3 .+ coeff5.*memslots5 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=2.9237775839655367 + 0.0im
    coeff2=0.0 + 0.04516040635859008im
    coeff3=0.1240818356655045 + 0.0im
    coeff4=0.0 - 0.019571570936427238im
    coeff5=2.425253007433925e-5 + 0.0im
    coeff6=1.0 + 0.0im
    # Smart lincomb recycle B5
    memslots7 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots3 .+ coeff5.*memslots5 .+ coeff6.*memslots7
    mul!(memslots7,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots6,memslots2,memslots7)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 7
    # Computing Ba7 = x*A+x*B2+x*B3+x*B4+x*B6
    coeff1=0.0 - 0.020637762737741035im
    coeff2=-1.093022784715649 + 0.0im
    coeff3=0.0 + 0.2537715581771087im
    coeff4=0.0005437426743473122 + 0.0im
    coeff5=1.0 + 0.0im
    memslots2 .= coeff1.*memslots1 .+ coeff2.*memslots4 .+ coeff3.*memslots3 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    # Computing Bb7 = x*A+x*B2+x*B3+x*B4+x*B6
    coeff1=0.0 - 0.020637762737741035im
    coeff2=-1.093022784715649 + 0.0im
    coeff3=0.0 + 0.2537715581771087im
    coeff4=0.0005437426743473122 + 0.0im
    coeff5=1.0 + 0.0im
    # Smart lincomb recycle A
    memslots1 .= coeff1.*memslots1 .+ coeff2.*memslots4 .+ coeff3.*memslots3 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 3
    # Deallocating B4 in slot 5
    # Deallocating B6 in slot 6
    # Computing B7 with operation: mult
    mul!(memslots3,memslots2,memslots1)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 1
    # Computing B8 with operation: mult
    mul!(memslots1,memslots3,memslots3)
    # Computing B9 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Deallocating B8 in slot 1
    # Computing B10 with operation: mult
    mul!(memslots1,memslots2,memslots2)
    # Deallocating B9 in slot 2
    # Computing B11 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Deallocating B10 in slot 1
    # Computing T2k13 = x*B7+x*B11
    coeff1=0.0 + 0.0im
    coeff2=1.0 + 0.0im
    # Smart lincomb recycle B7
    memslots3 .= coeff1.*memslots3 .+ coeff2.*memslots2
    # Deallocating B11 in slot 2
    return memslots3 # Returning T2k13
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
abstract type AbstractKoopmanAlgorithm end;

"""
    DMDPINV()

Approximates the koopman operator `K` based on

```julia
K = Y / X
```

where `Y` and `X` are data matrices.

"""
mutable struct DMDPINV <: AbstractKoopmanAlgorithm end;

# Fast but more allocations
(x::DMDPINV)(X::AbstractArray, Y::AbstractArray) = Y / X

"""
    DMDSVD(rtol)

Approximates the koopman operator `K` based on the singular value decomposition
of `X` such that

```julia
K = Y*V*Σ*U'
```

where `Y` and `X = U*Σ*V'` are data matrices.
If `rtol` ∈ (0, 1) is given, the singular value decomposition is reduced to include only
entries bigger than `rtol*maximum(Σ)`.
"""
mutable struct DMDSVD{T} <: AbstractKoopmanAlgorithm
    rtol::T
end;

DMDSVD() = DMDSVD(0.0)

# Slower but less allocations
function (x::DMDSVD)(X::AbstractArray, Y::AbstractArray)
    U, S, V = svd(X)
    if zero(x.rtol) < x.rtol < one(x.rtol)
        r = vec(S .> x.rtol*maximum(S))
        U = U[:, r]
        S = S[r]
        V = V[:, r]
    end

    return Y*V*Diagonal(one(eltype(X)) ./ S)*U'
end


text:
using SmartShockFinder
using Flux: onecold, crossentropy, @epochs, Data.DataLoader, params, Dropout, leakyrelu, softmax

datatyp = "NNPP"
n_dims = 2
η = 0.005                # learning rate
β = 0.005               # regularization paramater
number_epochs = 100
Sb = 500                 # batch size
L = 10                   # early stopping parameter
hidden_units_1 = 32
hidden_units_2 = 16
hidden_units_3 = 16
hidden_units_4 = 8
hidden_units_5 = 8

# traindata = joinpath(@__DIR__, "traindata2dNNPP.h5")
# download("https://gist.github.com/JuliaOd/f4d10153e539c026daf7076f087dc937/raw/traindata2dNNPP.h5", traindata)
# validdata = joinpath(@__DIR__, "validdata2dNNPP.h5")
# download("https://gist.github.com/JuliaOd/f4d10153e539c026daf7076f087dc937/raw/validdata2dNNPP.h5", validdata)

train_network2d(datatyp, n_dims, η, β, number_epochs, Sb, L,
                hidden_units_1, hidden_units_2, hidden_units_3,
                hidden_units_4, hidden_units_5)


text:
function fen_to_bitboard(fen_string::String)
    white = EMPTY
    R = zeros(UInt64, 0)
    N = zeros(UInt64, 0)
    B = zeros(UInt64, 0)
    Q = zeros(UInt64, 0)
    K = EMPTY
    P = zeros(UInt64, 0)

    black = EMPTY
    r = zeros(UInt64, 0)
    n = zeros(UInt64, 0)
    b = zeros(UInt64, 0)
    q = zeros(UInt64, 0)
    k = EMPTY
    p = zeros(UInt64, 0)

    free = EMPTY
    taken = EMPTY

    white_attacks = EMPTY
    black_attacks = EMPTY

    # board
    fen = split(fen_string, ' ')
    board = fen[1]
    square_i = 1
    fen_i = 1
    while true
        c = board[fen_i]
        if isnumeric(c)
            square_i += parse(Int64, c)
            fen_i += 1
        elseif c == '/'
            fen_i += 1
        else
            if isuppercase(c)
                # white
                white |= INT2UINT[square_i]
                if c == 'R'
                    push!(R, INT2UINT[square_i])
                elseif c == 'N'
                    push!(N, INT2UINT[square_i])
                elseif c == 'B'
                    push!(B, INT2UINT[square_i])
                elseif c == 'Q'
                    push!(Q, INT2UINT[square_i])
                elseif c == 'K'
                    K |= INT2UINT[square_i]
                elseif c == 'P'
                    push!(P, INT2UINT[square_i])
                end
            else
                # black
                black |= INT2UINT[square_i]
                if c == 'r'
                    push!(r, INT2UINT[square_i])
                elseif c == 'n'
                    push!(n, INT2UINT[square_i])
                elseif c == 'b'
                    push!(b, INT2UINT[square_i])
                elseif c == 'q'
                    push!(q, INT2UINT[square_i])
                elseif c == 'k'
                    k |= INT2UINT[square_i]
                elseif c == 'p'
                    push!(p, INT2UINT[square_i])
                end
            end
            square_i += 1
            fen_i += 1
        end

        if fen_i > length(board)

            if square_i != 65
                throw(ArgumentError("Invalid FEN string: too short/long"))
            end

            taken = white | black
            free = ~taken

            break
        end
    end

    if fen[2] == "w"
        player_color = "white"
    else
        player_color = "black"
    end

    white_king_moved = false
    black_king_moved = false
    if k != PGN2UINT["e8"]
        black_king_moved = true
    end
    if K != PGN2UINT["e1"]
        white_king_moved = true
    end
    white_can_castle_queenside = false
    white_can_castle_kingside = false
    black_can_castle_queenside = false
    black_can_castle_kingside = false
    for c in fen[3]
        if c == '-'
            break
        elseif c == 'K'
            white_can_castle_kingside = true
        elseif c == 'Q'
            white_can_castle_queenside = true
        elseif c == 'k'
            black_can_castle_kingside = true
        elseif c == 'q'
            black_can_castle_queenside = true
        end
    end

    enpassant_square = EMPTY
    enpassant_done = false
    if fen[4] != "-"
        enpassant_square = PGN2UINT[fen[4]]
    end

    if length(fen) >= 5
        halfmove_clock = parse(Int64, fen[5])
    else
        halfmove_clock = 1
        fullmove_clock = 0
    end
    if length(fen) == 6
        fullmove_clock = parse(Int64, fen[6])
    end

    game = Array{String,1}[]

    A = [EMPTY, EMPTY, EMPTY, EMPTY, EMPTY]
    a = [EMPTY, EMPTY, EMPTY, EMPTY, EMPTY]
    board = Bitboard(white, P, R, N, B, Q, K, A,
                    black, p, r, n, b, q, k, a,
                    free, taken,
                    white_attacks, black_attacks,
                    player_color,
                    white_can_castle_queenside,
                    white_can_castle_kingside,
                    black_can_castle_queenside,
                    black_can_castle_kingside,
                    white_king_moved, black_king_moved,
                    enpassant_square,
                    enpassant_done,
                    halfmove_clock, fullmove_clock,
                    fen_string, game)
    return update_attacked(board)
end


text:
module dpmeans
using Distributions, Random
using BSON
using Formatting: format

dropdim1(x) = dropdims(x, dims=1)
#= groupinds(x)
  for x isa Vector{Signed}.
  Extract the group of indices which match each unique value in x.
  Return, both the list of groups, and the unique elements they correspond to.
  adapted from GroupSlices.jl
 =#
function groupinds(ic::Vector{T}) where T <: Signed  # SignedInt
    d = Dict{T, T}()
    ia = unique(ic)
    n = length(ia)
    for i = 1:n
        d[ia[i]]= i
    end

    ib = Array{Vector{T}}(undef, n)
    for k = 1:n
        ib[k] = T[]
    end

    for h = 1:length(ic)
        push!(ib[d[ic[h]]], h)
    end
    return ib, ia
end

#=
  dist_to_centers.
  Helper function for dpmeans: calculate all distances ||x_i - μ_j|| over i, j.
 =#
function dist_to_centers(X::Matrix{T}, M::Matrix{T}; sq=false) where T <: Real
    n, d = size(X)
    k, _d2 = size(M)
    @assert (d == _d2) "X and M must have the same row dimension."

    X2    = sum(X .* X, dims=2)
    X_Mt  = X * M'
    M2    = sum(M .* M, dims=2)
    norm2 = (X2 .- 2X_Mt) .+ M2'
    if sq
        return norm2
    else
        return sqrt.(norm2)
    end
end

function dist_to_centers_colmajor(X::Matrix{T}, M::Matrix{T}; sq=false) where T <: Real
    d, n = size(X)
    _d2, k = size(M)
    @assert (d == _d2) "X and M must have the same row dimension."

    X2    = sum(X .* X, dims=1)
    X_Mt  = M' * X
    M2    = sum(M .* M, dims=1)
    norm2 = (X2 .- 2X_Mt) .+ M2'
    if sq
        return norm2
    else
        return sqrt.(norm2)
    end
end
#=
  row_mins.
  Helper function for dpmeans. Calculate the minimum of each row, and return
  the argmin in the second argument. MUCH MUCH faster than mapslices to
  `findmin`.
=#
function row_mins(X::Matrix{T}) where T <: Real
    n, d = size(X)
    out_v = Vector{T}(undef, n)
    out_ind = Vector{Int32}(undef, n)
    for nn = 1:n
        _min = X[nn, 1]
        _argmin = 1
        for dd = 2:d
            if X[nn, dd] < _min
                _min = X[nn, dd]
                _argmin = dd
            end
        end
        out_v[nn] = _min
        out_ind[nn] = _argmin
    end
    return out_v, out_ind
end

#=
  col_mins.
  Helper function for dpmeans. Calculate the minimum of each col, and return
  the argmin in the second argument. MUCH MUCH faster than mapslices to
  `findmin`.
=#
function col_mins(X::Matrix{T}) where T <: Real
    d, n = size(X)
    out_v = Vector{T}(undef, n)
    out_ind = Vector{Int32}(undef, n)
    for nn = 1:n
        _min = X[1, nn]
        _argmin = 1
        for dd = 2:d
            if X[dd, nn] < _min
                _min = X[dd, nn]
                _argmin = dd
            end
        end
        out_v[nn] = _min
        out_ind[nn] = _argmin
    end
    return out_v, out_ind
end

function pwise_mins(x::Vector{T}, y::Vector{T}, y_ind_id::Int32=Int32(2),
                    ind_if_x_wins::Union{Vector{Int32}, Nothing}=nothing,
                    val_if_x_wins::Union{Vector{T}, Nothing}=nothing) where T <: Real
    n = length(x)
    out_v = something(val_if_x_wins, copy(x))
    out_ind = something(ind_if_x_wins, ones(Int32, n))
    for nn = 1:n
        if x[nn] > y[nn]
            out_v[nn] = y[nn]
            out_ind[nn] = y_ind_id
        end
    end
    return out_v, out_ind
end

# function kpp_init(X::Array{Float64,2}, k::Int64)
#
#     n, d = size(X)
#     mu = zeros(k, d)
#     dist = Inf * ones(n)
#     λ = 1
#
#     idx = rand(1:n)
#     mu[1,:] = X[idx,:]
#     for i = 2:k
#         D = X .- mu[i-1,:]'
#         dist = min.(dist, sum(D.*D, dims=2))
#         idx_cond = rand() .< cumsum(dist/sum(dist))
#         idx_mu = findfirst(idx_cond)
#         mu[i,:] = X[idx_mu,:]
#         λ = maximum(dist)
#         println(λ)
#     end
#     return mu, λ
# end

# #= Batch version of dpmeans.
#     Nonparametric-like asymptotic σ->0 version of k-means with DP prior on k.
#     Loosely based / adapted from code originally written by Vadim Smolyakov.
#     https://github.com/vsmolyakov.
#     BATCH VERSION:
#     This actually works fairly poorly in general. A great counter example of why
#     one shouldn't do this is to imagine a dataset with a spherical distribution,
#     and initialising the first mean in the center. All points with distance λ from
#     the center will be moved to a new cluster, but the new cluster will also have
#     its mean.... at the center (or very close). Thus the routine will keep adding
#     new means without reducing the objective.
# =#
# function dpmeans_fit_batch(X::Array{T, 2}; k_init::Int64=3, max_iter::Int64=100) where T <: Number
#
#     #init params
#     k = k_init
#     n, d = size(X)
#
#     mu = mean(X, dims=1)
#     # display(mu)
#     # throw("Need to set lambda intelligently otherwise will continually add clusters.")
#     λ = 3
#     # mu, λ = kpp_init(X, k)
#
#     # println("lambda: ", λ)
#
#     obj = zeros(max_iter)
#     obj_tol = 1e-3
#
#     for iter = 1:max_iter
#         dist = zeros(n, k)
#
#         #assignment step
#         # custom functions are > order of magnitude faster
#         dist = dist_to_centers(X, mu; sq=true)
#         dmin, z = row_mins(dist)
#
#         #find(dist->dist==0, dist)
#         bad_idx = findall(dmin .> λ)
#         if length(bad_idx) > 0
#             k += 1
#             z[bad_idx] .= k
#             new_mean = mean(view(X, bad_idx, :), dims=1)
#             mu = vcat(mu, new_mean)
#             Xm = X .- new_mean
#             dist = hcat(dist, sum(Xm.*Xm, dims=2))
#         end
#
#         #update step
#         k_inds, ks = groupinds(z)
#         for j = 1:length(ks)
#             @views mu[ks[j],:] = mean(X[k_inds[j], :], dims=1)
#             @views obj[iter] += sum(dist[k_inds[j], ks[j]])
#         end
#         obj[iter] = obj[iter] + λ * k
#
#         #check convergence
#         if (iter > 1 && abs(obj[iter] - obj[iter-1]) < obj_tol * obj[iter])
#             println("converged in ", iter, " iterations.")
#             global Z = z
#             break
#         elseif iter == max_iter
#             @warn "DPmeans not converged"
#             global Z = z
#         end
#     end
#
#     return Z, mu, obj
# end

#= Gibbs version of dpmeans.
    Nonparametric-like asymptotic σ->0 version of k-means with DP prior on k.
    This version actually works, unlike the batch version! Recommend that the
    `shuffleX` option is used, as the algo is very sensitive to data order.
    Also there is a `collapse_thrsh` threshold which can remove small clusters
    if desired. I found these could be quite reasonable, and quite common, but
    annoying.

    > Expecting tall matrix: n * d. This is orders of magnitude slower than the
    default (wide) version for large matrices: it is not advised.
=#
function dpmeans_fit_tall(X::Matrix{T}; max_iter::Int=100, shuffleX::Bool=true,
                     lambda::T=1., collapse_thrsh::Int=0, suppress_warn=false) where T <: Number

    #init params
    k = 1
    n, d = size(X)
    @assert (n > d) "matrix must be arranged to be in tall format: n * d; found: tall format."

    # shuffling (and reordering for end)
    if shuffleX
        new_order = randperm(n)
        X = X[new_order, :]
        return_order(y::Vector) = y[sortperm(new_order)]
    else
        return_order = identity
    end

    mu = mean(X, dims=1)
    λ = lambda
    # mu, λ = kpp_init(X, k)

    obj = zeros(max_iter)
    nks = zeros(Int32, max_iter)

    obj_tol = 1e-3
    Z = Vector{Int}(undef, n)

    x_norm² = sum(X .* X, dims=2) |> _x -> dropdims(_x, dims=2)

    for iter = 1:max_iter

        # @timeit to "munorm" mu_norm² = sum(mu .* mu, dims=2)
        mu_terms = -2.0*X * mu' .+ sum(mu .* mu, dims=2)'
        c_mu_term_min, c_z = row_mins(mu_terms)

        for nn = 1:n
            dist = x_norm²[nn] + c_mu_term_min[nn]
            need_new_cls = dist > λ
            if need_new_cls
                k += 1
                # append new mean to mu
                mu = vcat(mu, X[nn,:]')
                addl_mu_terms = -2.0*X[nn:end,:] * X[nn,:] .+ sum(X[nn,:] .* X[nn,:])
                # update c_min / c_z
                c_mu_term_min[nn:end], c_z[nn:end] = pwise_mins(c_mu_term_min[nn:end], addl_mu_terms,
                                                                Int32(k), c_z[nn:end])
                dist = 0.
            end
            Z[nn] = c_z[nn]
            obj[iter] += dist
        end

        obj[iter] = obj[iter] + λ * k

        # update means
        k_inds, ks = groupinds(Z)
        for j = 1:length(ks)
            @views mu[ks[j],:] = mean(X[k_inds[j], :], dims=1)
        end
        nks[iter] = length(ks)

        #check convergence
        if (iter > 1 && abs(obj[iter] - obj[iter-1]) < obj_tol * obj[iter])
            # println("converged in ", iter, " iterations.")
            obj = obj[1:iter]
            nks = nks[1:iter]
            break
        elseif iter == max_iter
            !suppress_warn && @warn "DPmeans not converged"
        end
    end

    # collapse tiny clusters
    @info format("collapse thrsh is {:d}", collapse_thrsh)
    if collapse_thrsh > 0
        k_inds, ks = groupinds(Z)
        bad_ks_bool = map(length, k_inds) .< collapse_thrsh
        if sum(bad_ks_bool) > 0
            @info format("cluster(s) {:s} have <= 1 datapoint", findall(bad_ks_bool))
        end
        n_bad_k = sum(bad_ks_bool)
        if n_bad_k > 0
            bad_ixs = reduce(vcat, k_inds[findall(bad_ks_bool)])
            mu = mu[.!bad_ks_bool, :]

            dist = dist_to_centers(X[bad_ixs,:], mu; sq=true)
            _dmin, z = row_mins(dist)
            Z[bad_ixs] .= z
            nks[end] -= sum(bad_ks_bool)
            # not re-calculating obj because... who cares?
        end
    end
    return return_order(Z), mu, obj, nks
end

#= Gibbs version of dpmeans.
    Nonparametric-like asymptotic σ->0 version of k-means with DP prior on k.
    This version actually works, unlike the batch version! Recommend that the
    `shuffleX` option is used, as the algo is very sensitive to data order.
    Also there is a `collapse_thrsh` threshold which can remove small clusters
    if desired. I found these could be quite reasonable, and quite common, but
    annoying.
=#
function dpmeans_fit(X::Matrix{T}; max_iter::Int=100, shuffleX::Bool=true,
                     lambda::T=1., collapse_thrsh::Int=0, suppress_warn=false) where T <: Number

    #init params
    k = 1
    d, n = size(X)
    @assert (n > d) "matrix must be arranged to be in wide format: d * n; found: tall format."

    # shuffling (and reordering for end)
    if shuffleX
        new_order = randperm(n)
        X = X[:, new_order]
        return_order(y::Vector) = y[sortperm(new_order)]
    else
        return_order = identity
    end

    mu = mean(X, dims=2)
    λ = lambda

    obj = zeros(max_iter)
    nks = zeros(Int32, max_iter)

    obj_tol = 1e-3
    Z = Vector{Int}(undef, n)

    x_norm² = sum(X .* X, dims=1) |> dropdim1

    for iter = 1:max_iter

        # @timeit to "munorm" mu_norm² = sum(mu .* mu, dims=2)
        mu_terms = -2.0*mu'*X .+ sum(mu .* mu, dims=1)'
        c_mu_term_min, c_z = col_mins(mu_terms)

        for nn = 1:n
            dist = x_norm²[nn] + c_mu_term_min[nn]
            need_new_cls = dist > λ
            if need_new_cls
                k += 1
                # append new mean to mu
                mu = hcat(mu, X[:, nn])
                addl_mu_terms = -2.0*X[:, nn]' * X[:, nn:end] .+ sum(X[:, nn] .* X[:, nn]) |> dropdim1
                # update c_min / c_z
                c_mu_term_min[nn:end], c_z[nn:end] = pwise_mins(c_mu_term_min[nn:end], addl_mu_terms,
                                                                Int32(k), c_z[nn:end])
                dist = 0.
            end
            Z[nn] = c_z[nn]
            obj[iter] += dist
        end

        obj[iter] = obj[iter] + λ * k

        # update means
        k_inds, ks = groupinds(Z)
        for j = 1:length(ks)
            @views mu[:, ks[j]] = mean(X[:, k_inds[j]], dims=2)
        end
        nks[iter] = length(ks)

        #check convergence
        if (iter > 1 && abs(obj[iter] - obj[iter-1]) < obj_tol * obj[iter])
            # println("converged in ", iter, " iterations.")
            obj = obj[1:iter]
            nks = nks[1:iter]
            break
        elseif iter == max_iter
            !suppress_warn && @warn "DPmeans not converged"
        end
    end

    # collapse tiny clusters
    # bson("dbg/dpmean.bson", Z=Z, mu=mu, X=X)
    @info format("collapse thrsh is {:d}", collapse_thrsh)
    if collapse_thrsh > 0
        k_inds, ks = groupinds(Z)
        if minimum(ks) > 1   # only the first cluster can possibly be unused.
            mu = mu[:, 2:end]
            nks[end] -= 1
            k_inds, ks = groupinds(Z)
        end

        bad_ks_bool = map(length, k_inds) .< collapse_thrsh
        if sum(bad_ks_bool) > 0
            @info format("cluster(s) {:s} have <= 1 datapoint", findall(bad_ks_bool))
        end
        n_bad_k = sum(bad_ks_bool)
        # unlikely case that *all clusters* are too small
        !suppress_warn && n_bad_k == nks[end] && @warn "reducing collapse_thrsh as no clusters qualify"
        while n_bad_k == nks[end]
            collapse_thrsh *= 0.9
            bad_ks_bool = map(length, k_inds) .< collapse_thrsh
            n_bad_k = sum(bad_ks_bool)
        end
        if n_bad_k > 0
            bad_ixs = reduce(vcat, k_inds[findall(bad_ks_bool)])
            mu = mu[:, .!bad_ks_bool]

            dist = dist_to_centers_colmajor(X[:, bad_ixs], mu; sq=true)
            _dmin, z = col_mins(dist)
            Z[bad_ixs] .= z
            nks[end] -= sum(bad_ks_bool)
            # not re-calculating obj because... who cares?
        end
    end
    return return_order(Z), mu, obj, nks
end

end  # module


text:
struct SubMesh{N,T} <: AbstractMesh{N,T}
    mesh::GenericMesh{N,T}
    domain::Domain
    dom2elt::OrderedDict{DataType,Vector{Int}} # FIXME: this could be computed on the flight. Is that a problem?
    function SubMesh{N,T}(mesh::GenericMesh,Ω::Domain) where {N,T}
        idxs = dom2elt(mesh,Ω)
        return new{N,T}(mesh,Ω,idxs)
    end
end
SubMesh(m::GenericMesh{N,T},args...;kwargs...) where {N,T} = SubMesh{N,T}(m,args...;kwargs...)

Base.view(m::GenericMesh,Ω::Domain)           = SubMesh(m,Ω)
Base.view(m::GenericMesh,ent::AbstractEntity) = SubMesh(m,Domain(ent))

Base.getindex(m::GenericMesh,Ω::Domain) = view(m,Ω)
Base.getindex(m::GenericMesh,ent::AbstractEntity) = view(m,ent)

Base.getindex(m::SubMesh,Ω::Domain)           = view(parent(m),Ω)
Base.getindex(m::SubMesh,ent::AbstractEntity) = view(parent(m),ent)

parent(m::SubMesh) = m.mesh
domain(m::SubMesh) = m.domain

"""
    dom2elt(m::SubMesh,[E])

A dictionary with keys being the element types of `m`, and values being the
element indices in the parent mesh. If a type `E` is given, return the values
associated with that key.
"""
dom2elt(m::SubMesh) = m.dom2elt
dom2elt(m::SubMesh,E::Type{<:AbstractElement}) = m.dom2elt[E]

function etypes(submesh::SubMesh)
    Ω, M = submesh.domain, submesh.mesh
    ee = DataType[]
    for ent in entities(Ω)
        dict = M.ent2tags[ent]
        append!(ee, keys(dict))
    end
    return unique!(ee)
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
# ref: http://www.mit.edu/~yandros/doc/specs/fcgi-spec.html

BYTE(x,n) = UInt8((x>>((n-1)*8))&0xff)
B16(b1::UInt8, b0::UInt8) = (UInt16(b1) << 8) | b0

function padding(x::T) where {T<:Integer}
    rem = x % T(8)
    (rem > T(0)) ? (T(8) - rem) : T(0)
end

function _readnvlen(io::T) where {T <: IO}
    b = read(io, UInt8)
    len = UInt32(b)
    if (b >> 7) !== UInt8(0)
        len = ((len & 0x7f) << 24) | (UInt32(read(io, UInt8)) << 16) | (UInt32(read(io, UInt8)) << 8) | read(io, UInt8)
    end
    len
end

function _writenvlen(io::I, l::T) where {I<:IO,T<:Integer}
    nbytes = 0
    if l > 0x7f
        l32 = UInt32(l)
        nbytes += write(io, BYTE(l32, 4) | 0x80)
        nbytes += write(io, BYTE(l32, 3))
        nbytes += write(io, BYTE(l32, 2))
        nbytes += write(io, BYTE(l32, 1))
    else
        nbytes += write(io, UInt8(l))
    end
    nbytes
end

#=
function _readnvlen(d::Vector{UInt8}, pos::Int)
    len = UInt32(d[pos])
    pos += 1
    if (len >> 7) !== 0
        len = ((len & 0x7f) << 24)
        pos += 1
        len |= (UInt32(d[pos]) << 16)
        pos += 1
        len |= (UInt32(d[pos]) << 8)
        pos += 1
        len |= UInt32(d[pos])
    end
    len, pos
end
=#

function _readbytes(io::I, len::T) where {I<:IO, T<:Integer}
    @assert len >= 0
    (len == 0) && (return UInt8[])
    buff = read(io, len)
    while length(buff) != len
        append!(buff, read(io, (len - length(buff))))
    end
    buff
end
_readstr(io::I, len::T) where {I<:IO, T<:Integer} = String(_readbytes(io, len))

"""Listening socket file number"""
const LISTENSOCK_FILENO = 0
const MAX_LENGTH = 0xffff

"""Number of bytes in a FCGIHeader.  Future versions of the protocol will not reduce this number."""
const HEADER_LEN = UInt8(8)

"""Value for version component of FCGIHeader"""
const VERSION_1 = UInt8(1)

"""Value for requestId component of FCGIHeader"""
const NULL_REQUEST_ID = 0

"""Values for type component of FCGIHeader"""
struct _FCGIHeaderType
    BEGIN_REQUEST::UInt8        # FCGIBeginRequest
    ABORT_REQUEST::UInt8        # not applicable
    END_REQUEST::UInt8          # FCGIEndRequest
    PARAMS::UInt8               # FCGIParams
    STDIN::UInt8                # FCGIRecord.content
    STDOUT::UInt8               # FCGIRecord.content
    STDERR::UInt8               # FCGIRecord.content
    DATA::UInt8                 # FCGIRecord.content (applicable only in context of filters)
    GET_VALUES::UInt8           # FCGIParams with empty values
    GET_VALUES_RESULT::UInt8    # FCGIParams
    UNKNOWN_TYPE::UInt8         # FCGIUnknownType
    MAXTYPE::UInt8
    function _FCGIHeaderType()
        new(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 11)
    end
end
const FCGIHeaderType = _FCGIHeaderType()
function reqtypetostring(type::UInt8)
    (type === FCGIHeaderType.BEGIN_REQUEST) ? "BEGIN_REQUEST" :
    (type === FCGIHeaderType.ABORT_REQUEST) ? "ABORT_REQUEST" :
    (type === FCGIHeaderType.END_REQUEST) ? "END_REQUEST" :
    (type === FCGIHeaderType.PARAMS) ? "PARAMS" :
    (type === FCGIHeaderType.STDIN) ? "STDIN" :
    (type === FCGIHeaderType.STDOUT) ? "STDOUT" :
    (type === FCGIHeaderType.STDERR) ? "STDERR" :
    (type === FCGIHeaderType.DATA) ? "DATA" :
    (type === FCGIHeaderType.GET_VALUES) ? "GET_VALUES" :
    (type === FCGIHeaderType.GET_VALUES_RESULT) ? "GET_VALUES_RESULT" : "UNKNOWN"
end

struct FCGIHeader
    version::UInt8
    type::UInt8
    requestIdB1::UInt8
    requestIdB0::UInt8
    contentLengthB1::UInt8
    contentLengthB0::UInt8
    paddingLength::UInt8
    # reserved::UInt8

    FCGIHeader(io::T) where {T <: IO} = FCGIHeader(read(io, UInt64))
    FCGIHeader(b::UInt64) = new(BYTE(b,1), BYTE(b,2), BYTE(b,3), BYTE(b,4), BYTE(b,5), BYTE(b,6), BYTE(b,7))
    function FCGIHeader(type::UInt8, reqid::UInt16, contentlen::UInt16)
        new(VERSION_1, type, BYTE(reqid,2), BYTE(reqid,1), BYTE(contentlen,2), BYTE(contentlen,1), padding(contentlen))
    end
end
requestid(hdr::FCGIHeader) = B16(hdr.requestIdB1, hdr.requestIdB0)
contentlength(hdr::FCGIHeader) = B16(hdr.contentLengthB1, hdr.contentLengthB0)
show(io::IO, hdr::FCGIHeader) = print(io, "FCGIHeader(", reqtypetostring(hdr.type), ", id=", requestid(hdr), ", len=", contentlength(hdr), ")")

"""Mask for flags component of FCGIBeginRequest"""
const KEEP_CONN = UInt8(1)

"""Values for role component of FCGIBeginRequest"""
struct _FCGIRequestRole
    RESPONDER::UInt16
    AUTHORIZER::UInt16
    FILTER::UInt16
    function _FCGIRequestRole()
        new(1, 2, 3)
    end
end
const FCGIRequestRole = _FCGIRequestRole()

struct FCGIBeginRequest
    roleB1::UInt8
    roleB0::UInt8
    flags::UInt8
    # reserved::Vector{UInt8} # length 5
    FCGIBeginRequest(d::Vector{UInt8}) = new(d[1], d[2], d[3])
    function FCGIBeginRequest(role::UInt16, flags::UInt8=UInt8(0))
        new(BYTE(role,2), BYTE(role,1), flags)
    end
end
role(req::FCGIBeginRequest) = B16(req.roleB1, req.roleB0)
keepconn(req::FCGIBeginRequest) = (req.flags & KEEP_CONN) != 0x00

"""Values for protocolStatus component of FCGIEndRequest"""
struct _FCGIEndRequestProtocolStatus
    REQUEST_COMPLETE::UInt8
    CANT_MPX_CONN::UInt8
    OVERLOADED::UInt8
    UNKNOWN_ROLE::UInt8
    function _FCGIEndRequestProtocolStatus()
        new(0, 1, 2, 3)
    end
end
const FCGIEndRequestProtocolStatus = _FCGIEndRequestProtocolStatus()

struct FCGIEndRequest
    appStatusB3::UInt8
    appStatusB2::UInt8
    appStatusB1::UInt8
    appStatusB0::UInt8
    protocolStatus::UInt8
    # reserved::Vector{UInt8} # length 3
    FCGIEndRequest(d::Vector{UInt8}) = new(d[1], d[2], d[3], d[4], d[5])
    function FCGIEndRequest(appstatus::UInt32, protocolstatus::UInt8=FCGIEndRequestProtocolStatus.REQUEST_COMPLETE)
        new(BYTE(appstatus,4), BYTE(appstatus,3), BYTE(appstatus,2), BYTE(appstatus,1), protocolstatus)
    end
end

# Variable names for FCGI_GET_VALUES / FCGI_GET_VALUES_RESULT records
"""The maximum number of concurrent transport connections this application will accept, e.g. "1" or "10"."""
const MAX_CONNS  = "FCGI_MAX_CONNS"
"""The maximum number of concurrent requests this application will accept, e.g. "1" or "50"."""
const MAX_REQS   = "FCGI_MAX_REQS"
"""If this application does not multiplex connections (i.e. handle concurrent requests over each connection) then "0", "1" otherwise."""
const MPXS_CONNS = "FCGI_MPXS_CONNS"

struct FCGIUnknownType
    type::UInt8
    # reserved::Vector{UInt8} # length 7
end
FCGIUnknownType(d::Vector{UInt8}) = FCGIUnknownType(d[1])

struct FCGINameValuePair
    name::String
    value::String
end
function FCGINameValuePair(io::T) where {T <: IO}
    lname = _readnvlen(io)
    lvalue = _readnvlen(io)
    name = _readstr(io, lname)
    value = _readstr(io, lvalue)
    FCGINameValuePair(name, value)
end

struct FCGIParams
    nvpairs::Vector{FCGINameValuePair}

    FCGIParams() = new(Vector{FCGINameValuePair}())
    FCGIParams(d::Vector{UInt8}) = FCGIParams(IOBuffer(d))
    function FCGIParams(io::T) where {T <: IO}
        nvpairs = Vector{FCGINameValuePair}()
        while !eof(io)
            push!(nvpairs, FCGINameValuePair(io))
        end
        new(nvpairs)
    end
    function FCGIParams(nv::Dict{String,String})
        nvpairs = [FCGINameValuePair(n,v) for (n,v) in nv]
        new(nvpairs)
    end
end

struct FCGIRecord
    header::FCGIHeader
    content::Vector{UInt8}

    function FCGIRecord(io::T) where {T <: IO}
        header = FCGIHeader(io)
        content = _readbytes(io, contentlength(header))
        _ = _readbytes(io, header.paddingLength) # discard padding
        new(header, content)
    end
    function FCGIRecord(type::UInt8, reqid::UInt16, content::Vector{UInt8})
        contentlen = UInt16(length(content))
        new(FCGIHeader(type, reqid, contentlen), content)
    end
    FCGIRecord(type::UInt8, reqid::UInt16, content::String) = FCGIRecord(type, reqid, convert(Vector{UInt8},codeunits(content)))
    function FCGIRecord(type::UInt8, reqid::UInt16, content::Union{FCGIBeginRequest,FCGIEndRequest,FCGIParams,FCGIUnknownType})
        iob = IOBuffer()
        fcgiwrite(iob, content)
        FCGIRecord(type, reqid, take!(iob))
    end
end
show(io::IO, rec::FCGIRecord) = print(io, "FCGIRecord(", reqtypetostring(rec.header.type), ", id=", requestid(rec.header), ", len=", contentlength(rec.header), ")")

fcgiwrite(io::T, hdr::FCGIHeader) where {T<:IO} = write(io, hdr.version, hdr.type, hdr.requestIdB1, hdr.requestIdB0, hdr.contentLengthB1, hdr.contentLengthB0, hdr.paddingLength, UInt8(0))
fcgiwrite(io::T, begreq::FCGIBeginRequest) where {T<:IO} = write(io, begreq.roleB1, begreq.roleB0, begreq.flags, zeros(UInt8, 5))
fcgiwrite(io::T, endreq::FCGIEndRequest) where {T<:IO} = write(io, endreq.appStatusB3, endreq.appStatusB2, endreq.appStatusB1, endreq.appStatusB0, endreq.protocolStatus, zeros(UInt8, 3))
fcgiwrite(io::T, unk::FCGIUnknownType) where {T<:IO} = write(io, unk.type, zeros(UInt8, 7))
function fcgiwrite(io::T, nv::FCGINameValuePair) where {T<:IO}
    nbytes = 0
    lname = length(nv.name)
    lvalue = length(nv.value)
    nbytes += _writenvlen(io, lname)
    nbytes += _writenvlen(io, lvalue)
    nbytes += write(io, nv.name, nv.value)
    nbytes
end
fcgiwrite(io::T, params::FCGIParams) where {T<:IO} = sum([fcgiwrite(io, nv) for nv in params.nvpairs])
function fcgiwrite(io::T, rec::FCGIRecord) where {T<:IO}
    iob = IOBuffer()
    nbytes = fcgiwrite(iob, rec.header)
    nbytes += write(iob, rec.content)
    if rec.header.paddingLength > 0
        nbytes += write(iob, zeros(UInt8, rec.header.paddingLength))
    end
    writefully(io, iob)
end

const FCGIMsgComponent = Union{FCGIBeginRequest, FCGIEndRequest, FCGIUnknownType, FCGINameValuePair, FCGIParams, FCGIRecord}
function ==(o1::T, o2::T) where {T <: FCGIMsgComponent}
    for fld in fieldnames(T)
        (getfield(o1, fld) == getfield(o2, fld)) || (return false)
    end
    true
end

const FCGIServerSocket = Union{Sockets.TCPServer,Sockets.PipeServer}
const FCGISocket = Union{TCPSocket,Base.PipeEndpoint}


text:
#=
setindex.jl

Provide setindex!() capabilities like A[i,j] = s for LinearMapAM objects.

The setindex! function here is (necessarily) quite inefficient.
It is provided solely so that a LinearMapAM conforms to the AbstractMatrix
requirement that a setindex! method exists.
Use of this function is strongly discouraged, other than for testing.

2018-01-19, Jeff Fessler, University of Michigan
=#

using LinearMaps #: LinearMaps.FunctionMap

#Indexer = AbstractVector{Int}

"""
`A[i,j] = X`

Mathematically, if B = copy(A) and then we say `B[i,j] = s`, for scalar `s`,
then `B = A + (s - A[i,j]) e_i e_j'`
where `e_i` and `e_j` are the appropriate unit vectors.
Using this basic math, we can perform `B*x` using `A*x` and `A[i,j]`.

This idea generalizes to `B[ii,jj] = X` where `ii` and `jj` are vectors.

This method works because LinearMapAM is a *mutable* struct.
"""
function Base.setindex!(A::LinearMapAM, X::AbstractMatrix,
    ii::Indexer, jj::Indexer,
)

    # todo: handle WrappedMap differently
    L = A._lmap
    Aij = A[ii,jj] # must extract this outside of forw() to avoid recursion!
    Aij = reshape(Aij, size(X)) # needed for [i,:] case
    if true
        forw = (x) ->
            begin
                tmp = L * x
                tmp[ii] += (X - Aij) * x[jj]
                return tmp
            end
    end

    has_adj = !(typeof(L) <: LinearMaps.FunctionMap) || (L.fc !== nothing)
    if has_adj
        back = (y) ->
            begin
                tmp = L' * y
                tmp[jj] += (X - Aij)' * y[ii]
                return tmp
            end
    end

#=
    # for future reference, some of this could be implemented
    A.issymmetric = false # could check if i==j
    A.ishermitian = false # could check if i==j and v=v'
    A.posdef = false # could check if i==j and v > A[i,j]
=#

    if has_adj
        A._lmap = LinearMap(forw, back, size(L)...)
    else
        A._lmap = LinearMap(forw, size(L)...)
    end

    nothing
end

# [i,j] = s
Base.setindex!(A::LinearMapAM, s::Number, i::Int, j::Int) =
    setindex!(A, fill(s,1,1), [i], [j])

# [ii,jj] = X
Base.setindex!(A::LinearMapAM, X::AbstractMatrix,
    ii::AbstractVector{Bool}, jj::AbstractVector{Bool}) =
    setindex!(A, X, findall(ii), findall(jj))

# [:,jj] = X
Base.setindex!(A::LinearMapAM, X::AbstractMatrix, ::Colon, jj::AbstractVector{Bool}) =
    setindex!(A, X, :, findall(jj))
Base.setindex!(A::LinearMapAM, X::AbstractMatrix, ::Colon, jj::Indexer) =
    setindex!(A, X, 1:size(A,1), jj)

# [ii,:] = X
Base.setindex!(A::LinearMapAM, X::AbstractMatrix, ii::AbstractVector{Bool}, ::Colon) =
    setindex!(A, X, findall(ii), :)
Base.setindex!(A::LinearMapAM, X::AbstractMatrix, ii::Indexer, ::Colon) =
    setindex!(A, X, ii, 1:size(A,2))

# [:,j] = v
Base.setindex!(A::LinearMapAM, v::AbstractVector, ::Colon, j::Int) =
    setindex!(A, reshape(v,:,1), :, [j])

# [ii,j] = v
Base.setindex!(A::LinearMapAM, v::AbstractVector, ii::Indexer, j::Int) =
    setindex!(A, reshape(v,:,1), ii, [j])

# [i,:] = v
Base.setindex!(A::LinearMapAM, v::AbstractVector, i::Int, ::Colon) =
    setindex!(A, reshape(v,1,:), [i], :)

# [i,jj] = v
Base.setindex!(A::LinearMapAM, v::AbstractVector, i::Int, jj::Indexer) =
    setindex!(A, reshape(v,1,:), [i], jj)

# [ii,jj] = s
Base.setindex!(A::LinearMapAM, s::Number, ii::AbstractVector{Bool}, jj::Indexer) =
    setindex!(A, s, findall(ii), jj)
Base.setindex!(A::LinearMapAM, s::Number, ii::Indexer, jj::AbstractVector{Bool}) =
    setindex!(A, s, ii, findall(jj))
Base.setindex!(A::LinearMapAM, s::Number,
    ii::AbstractVector{Bool}, jj::AbstractVector{Bool}) =
    setindex!(A, s, findall(ii), findall(jj))
Base.setindex!(A::LinearMapAM, s::Number, ii::Indexer, jj::Indexer) =
    setindex!(A, fill(s,length(ii),length(jj)), ii, jj)

# [:,j] = s
Base.setindex!(A::LinearMapAM, s::Number, ::Colon, j::Int) =
    setindex!(A, s, 1:size(A,1), [j])

# [ii,:] = s
Base.setindex!(A::LinearMapAM, s::Number, ii::AbstractVector{Bool}, ::Colon) =
    setindex!(A, s, findall(ii), :)
Base.setindex!(A::LinearMapAM, s::Number, ii::Indexer, ::Colon) =
    setindex!(A, s, ii, 1:size(A,2))

# [ii,j] = s
Base.setindex!(A::LinearMapAM, s::Number, ii::AbstractVector{Bool}, j::Int) =
    setindex!(A, s, findall(ii), [j])
Base.setindex!(A::LinearMapAM, s::Number, ii::Indexer, j::Int) =
    setindex!(A, fill(s, length(ii), 1), ii, [j])

# [i,:] = s
Base.setindex!(A::LinearMapAM, s::Number, i::Int, ::Colon) =
    setindex!(A, fill(s, 1, size(A,2)), [i], :)

# [i,jj] = s
Base.setindex!(A::LinearMapAM, s::Number, i::Int, jj::AbstractVector{Bool}) =
    setindex!(A, s, i, findall(jj))
Base.setindex!(A::LinearMapAM, s::Number, i::Int, jj::Indexer) =
    setindex!(A, fill(s, 1, length(jj)), [i], jj)

# [:,jj] = s
Base.setindex!(A::LinearMapAM, s::Number, ::Colon, jj::AbstractVector{Bool}) =
    setindex!(A, s, :, findall(jj))
Base.setindex!(A::LinearMapAM, s::Number, ::Colon, jj::Indexer) =
    setindex!(A, s, 1:size(A,1), jj)

# [kk]
#= too much work, so unsupported
Base.setindex!(A::LinearMapAM, v::AbstractVector,
        kk::Indexer) =
=#

# [:] = v
Base.setindex!(A::LinearMapAM, v::AbstractVector, ::Colon) =
    setindex!(A, reshape(v,size(A)), :, :)

# [:] = s
function Base.setindex!(A::LinearMapAM, s::Number, ::Colon)
    (M,N) = size(A)
    forw = x -> fill(s, M) * (ones(1,N) * x)[1]
    back = y -> fill(conj(s), N) * (ones(1, M) * y)[1]
    A._lmap = LinearMap(forw, back, M, N)
end

# [:,:] = s
Base.setindex!(A::LinearMapAM, s::Number, ::Colon, ::Colon) =
    setindex!(A, s, :)

# [:,:] = X
function Base.setindex!(A::LinearMapAM, X::AbstractMatrix, ::Colon, ::Colon)
    A._lmap = LinearMap(X)
end

# [k] = s
"""
`setindex!(A::LinearMapAM, s::Number, k::Int)`

Provide the single index version to meet the `AbstractArray` spec:
https://docs.julialang.org/en/latest/manual/interfaces/#Indexing-1
"""
Base.setindex!(A::LinearMapAM, s::Number, k::Int) =
    setindex!(A, s, Tuple(CartesianIndices(size(A))[k])...) # any better way?

#=
"""
`setindex!(A::LinearMapAM, v::Number, k::Int)`

Provide the single index version to meet the `AbstractArray` spec:
https://docs.julialang.org/en/latest/manual/interfaces/#Indexing-1
"""
function Base.setindex!(A::LinearMapAM, v::Number, k::Int)
    c = CartesianIndices(size(A))[k] # is there a more elegant way?
    setindex!(A, v::Number, c[1], c[2])
end
=#


text:
using BinaryBuilder

cuda_version = v"11.1.0"   # NOTE: could be less specific

sources_linux_x64 = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.4/11.1_20200923/cudnn-11.1-linux-x64-v8.0.4.30.tgz",
                  "8f4c662343afce5998ce963500fe3bb167e9a508c1a1a949d821a4b80fa9beab")
]
sources_linux_ppc64le = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.4/11.1_20200923/cudnn-11.1-linux-ppc64le-v8.0.4.30.tgz",
                  "b4ddb51610cbae806017616698635a9914c3e1eb14259f3a39ee5c84e7106712")
]
sources_windows = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.4/11.1_20200923/cudnn-11.1-windows-x64-v8.0.4.30.zip",
                  "58bf1fb324d11088a28c3a9f71a38f5dfe167efdf723815dce1867bc03ddaea2")
]

include("../common.jl")