text:
using Test, LinearAlgebra
using DiffEqSensitivity, StochasticDiffEq
using ForwardDiff, Zygote
using Random

@info "SDE Non-Diagonal Noise Adjoints"

seed = 100
Random.seed!(seed)

tstart = 0.0
tend = 0.1
dt = 0.005
trange = (tstart, tend)
t = tstart:dt:tend
tarray = collect(t)

function g(u,p,t)
  sum(u.^2.0/2.0)
end

function dg!(out,u,p,t,i)
  (out.=-u)
end

# non-diagonal noise
@testset "Non-diagonal noise tests" begin
    Random.seed!(seed)

    u₀ = [0.75,0.5]
    p = [-1.5,0.05,0.2, 0.01]

    dtnd = tend/1e3

    # Example from Roessler, SIAM J. NUMER. ANAL, 48, 922–952 with d = 2; m = 2
    function f_nondiag!(du,u,p,t)
      du[1] = p[1]*u[1] + p[2]*u[2]
      du[2] = p[2]*u[1] + p[1]*u[2]
      nothing
    end

    function g_nondiag!(du,u,p,t)
      du[1,1] = p[3]*u[1] + p[4]*u[2]
      du[1,2] = p[3]*u[1] + p[4]*u[2]
      du[2,1] = p[4]*u[1] + p[3]*u[2]
      du[2,2] = p[4]*u[1] + p[3]*u[2]
      nothing
    end

    function f_nondiag(u,p,t)
      dx = p[1]*u[1] + p[2]*u[2]
      dy = p[2]*u[1] + p[1]*u[2]
      [dx,dy]
    end

    function g_nondiag(u,p,t)
      du11 = p[3]*u[1] + p[4]*u[2]
      du12 = p[3]*u[1] + p[4]*u[2]
      du21 = p[4]*u[1] + p[3]*u[2]
      du22 = p[4]*u[1] + p[3]*u[2]

      [du11  du12
        du21  du22]
    end

    function f_nondiag_analytic(u0,p,t,W)
      A = [[p[1], p[2]] [p[2], p[1]]]
      B = [[p[3], p[4]] [p[4], p[3]]]
      tmp = A*t + B*W[1] + B*W[2]
      exp(tmp)*u0
    end

    noise_matrix = similar(p,2,2)
    noise_matrix .= false

    Random.seed!(seed)
    prob = SDEProblem(f_nondiag!,g_nondiag!,u₀,trange,p,noise_rate_prototype=noise_matrix)
    sol = solve(prob, EulerHeun(), dt=dtnd, save_noise=true)

    noise_matrix = similar(p,2,2)
    noise_matrix .= false
    Random.seed!(seed)
    proboop = SDEProblem(f_nondiag,g_nondiag,u₀,trange,p,noise_rate_prototype=noise_matrix)
    soloop = solve(proboop,EulerHeun(), dt=dtnd, save_noise=true)

    res_sde_u0, res_sde_p = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=BacksolveAdjoint())

    @info res_sde_p

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ReverseDiffNoise()))

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=BacksolveAdjoint(noise=false))

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ZygoteNoise()))

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-4)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ReverseDiffNoise()))

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-4)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=InterpolatingAdjoint(noise=false))

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-4)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=BacksolveAdjoint())

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ZygoteNoise()))

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=BacksolveAdjoint(noise=false))

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=InterpolatingAdjoint())

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-4)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ZygoteNoise()))

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-4)

    @info res_sde_pa

    res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
        ,dt=dtnd,adaptive=false,sensealg=InterpolatingAdjoint(noise=false))

    @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
    @test isapprox(res_sde_pa, res_sde_p, rtol=1e-4)

    @info res_sde_pa

    function compute_grads_nd(sol)
      xdis = sol(tarray)

      mat1 = Matrix{Int}(I, 2, 2)
      mat2 = ones(2,2)-mat1

      tmp1 = similar(p)
      tmp1 *= false

      tmp2 = similar(xdis.u[1])
      tmp2 *= false

      for (i, u) in enumerate(xdis)
        tmp1[1]+=xdis.t[i]*u'*mat1*u
        tmp1[2]+=xdis.t[i]*u'*mat2*u
        tmp1[3]+=sum(sol.W(xdis.t[i])[1])*u'*mat1*u
        tmp1[4]+=sum(sol.W(xdis.t[i])[1])*u'*mat2*u

        tmp2 += u.^2
      end

      return tmp2 ./ xdis.u[1], tmp1
    end

    res1, res2 = compute_grads_nd(soloop)

    @test isapprox(res1, res_sde_u0, rtol=1e-4)
    @test isapprox(res2, res_sde_p', rtol=1e-4)

end

@testset "diagonal but mixing noise tests" begin

  Random.seed!(seed)
  u₀ = [0.75,0.5]
  p = [-1.5,0.05,0.2, 0.01]
  dtmix = tend/1e3

  # Example from Roessler, SIAM J. NUMER. ANAL, 48, 922–952 with d = 2; m = 2
  function f_mixing!(du,u,p,t)
    du[1] = p[1]*u[1] + p[2]*u[2]
    du[2] = p[2]*u[1] + p[1]*u[2]
    nothing
  end

  function g_mixing!(du,u,p,t)
    du[1] = p[3]*u[1] + p[4]*u[2]
    du[2] = p[3]*u[1] + p[4]*u[2]
    nothing
  end

  function f_mixing(u,p,t)
    dx = p[1]*u[1] + p[2]*u[2]
    dy = p[2]*u[1] + p[1]*u[2]
    [dx,dy]
  end

  function g_mixing(u,p,t)
    dx = p[3]*u[1] + p[4]*u[2]
    dy = p[3]*u[1] + p[4]*u[2]
    [dx,dy]
  end

  Random.seed!(seed)
  prob = SDEProblem(f_mixing!,g_mixing!,u₀,trange,p)

  soltsave = collect(trange[1]:dtmix:trange[2])
  sol = solve(prob, EulerHeun(), dt=dtmix, save_noise=true, saveat=soltsave)

  Random.seed!(seed)
  proboop = SDEProblem(f_mixing,g_mixing,u₀,trange,p)
  soloop = solve(proboop,EulerHeun(), dt=dtmix, save_noise=true, saveat=soltsave)

  #oop

  res_sde_u0, res_sde_p = adjoint_sensitivities(soloop,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noisemixing=true))

  @info res_sde_p

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
    ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ZygoteNoise(), noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
    ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noise=false,noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,Array(t)
      ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ReverseDiffNoise(), noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noisemixing=true, noise=DiffEqSensitivity.ZygoteNoise()))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noise=false, noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5)

  @info res_sde_pa

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(soloop,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ReverseDiffNoise(),noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5)

  @info res_sde_pa

  res_sde_u0, res_sde_p = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
      ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5)

  @info res_sde_p

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
      ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ZygoteNoise(), noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
      ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noise=false,noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
      ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ReverseDiffNoise(), noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-6)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-6)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
      ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5) # would pass with 1e-4 but last noise value is off

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
      ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noisemixing=true, noise=DiffEqSensitivity.ZygoteNoise()))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noise=false, noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5)

  @info res_sde_pa

  res_sde_u0a, res_sde_pa = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
      ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ReverseDiffNoise(),noisemixing=true))

  @test isapprox(res_sde_u0a, res_sde_u0, rtol=1e-5)
  @test isapprox(res_sde_pa, res_sde_p, rtol=1e-5)

  @info res_sde_pa

  function GSDE(p)
    Random.seed!(seed)
    tmp_prob = remake(prob,u0=eltype(p).(prob.u0),p=p,
                  tspan=eltype(p).(prob.tspan))
    _sol = solve(tmp_prob,EulerHeun(),dt=dtmix,adaptive=false,saveat=Array(t))
    A = convert(Array,_sol)
    res = g(A,p,nothing)
  end

  res_sde_forward = ForwardDiff.gradient(GSDE,p)

  @test isapprox(res_sde_p', res_sde_forward, rtol=1e-5)

  function GSDE2(u0)
    Random.seed!(seed)
    tmp_prob = remake(prob,u0=u0,p=eltype(p).(prob.p),
                      tspan=eltype(p).(prob.tspan))
    _sol = solve(tmp_prob,EulerHeun(),dt=dtmix,adaptive=false,saveat=Array(t))
    A = convert(Array,_sol)
    res = g(A,p,nothing)
  end

  res_sde_forward = ForwardDiff.gradient(GSDE2,u₀)

  @test isapprox(res_sde_forward, res_sde_u0, rtol=1e-5)

end

@testset "mixing noise inplace/oop tests" begin
  Random.seed!(seed)
  u₀ = [0.75,0.5]
  p = [-1.5,0.05,0.2, 0.01]
  dtmix = tend/1e3

  # Example from Roessler, SIAM J. NUMER. ANAL, 48, 922–952 with d = 2; m = 2
  function f_mixing!(du,u,p,t)
    du[1] = p[1]*u[1] + p[2]*u[2]
    du[2] = p[2]*u[1] + p[1]*u[2]
    nothing
  end

  function g_mixing!(du,u,p,t)
    du[1] = p[3]*u[1] + p[4]*u[2]
    du[2] = p[3]*u[1] + p[4]*u[2]
    nothing
  end

  function f_mixing(u,p,t)
    dx = p[1]*u[1] + p[2]*u[2]
    dy = p[2]*u[1] + p[1]*u[2]
    [dx,dy]
  end

  function g_mixing(u,p,t)
    dx = p[3]*u[1] + p[4]*u[2]
    dy = p[3]*u[1] + p[4]*u[2]
    [dx,dy]
  end

  Random.seed!(seed)
  prob = SDEProblem(f_mixing!,g_mixing!,u₀,trange,p)
  soltsave = collect(trange[1]:dtmix:trange[2])
  sol = solve(prob, EulerHeun(), dt=dtmix, save_noise=true, saveat=soltsave)

  Random.seed!(seed)
  proboop = SDEProblem(f_mixing,g_mixing,u₀,trange,p)
  soloop = solve(proboop,EulerHeun(), dt=dtmix, save_noise=true, saveat=soltsave)

  @test sol.u ≈ soloop.u atol = 1e-14

  # BacksolveAdjoint

  res_sde_u0, res_sde_p = adjoint_sensitivities(soloop,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noisemixing=true))

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=BacksolveAdjoint(noisemixing=true))

  @test res_sde_u0 ≈ res_sde_u02 atol = 1e-14
  @test res_sde_p ≈ res_sde_p2 atol = 1e-14

  @show res_sde_u0

  adjproboop = SDEAdjointProblem(soloop,BacksolveAdjoint(noisemixing=true),dg!,tarray, nothing)
  adj_soloop = solve(adjproboop,EulerHeun(); dt=dtmix, tstops=soloop.t, adaptive=false)

  @test adj_soloop[end][length(p)+length(u₀)+1:end] == soloop.u[1]
  @test - adj_soloop[end][1:length(u₀)] == res_sde_u0
  @test adj_soloop[end][length(u₀)+1:end-length(u₀)] == res_sde_p'

  adjprob = SDEAdjointProblem(sol,BacksolveAdjoint(noisemixing=true, checkpointing=true),dg!,tarray, nothing)
  adj_sol = solve(adjprob,EulerHeun(); dt=dtmix, adaptive=false,tstops=soloop.t)

  @test adj_soloop[end] ≈  adj_sol[end]  rtol=1e-15

  adjprob = SDEAdjointProblem(sol,BacksolveAdjoint(noisemixing=true, checkpointing=false),dg!,tarray, nothing)
  adj_sol = solve(adjprob,EulerHeun(); dt=dtmix, adaptive=false,tstops=soloop.t)

  @test adj_soloop[end] ≈  adj_sol[end]  rtol=1e-8

  # InterpolatingAdjoint

  res_sde_u0, res_sde_p = adjoint_sensitivities(soloop,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noisemixing=true))

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol,EulerHeun(),dg!,tarray
    ,dt=dtmix,adaptive=false,sensealg=InterpolatingAdjoint(noisemixing=true))

  @test res_sde_u0 ≈ res_sde_u02 atol = 1e-8
  @test res_sde_p ≈ res_sde_p2 atol = 5e-8

  @show res_sde_u0

  adjproboop = SDEAdjointProblem(soloop,InterpolatingAdjoint(noisemixing=true),dg!,tarray, nothing)
  adj_soloop = solve(adjproboop,EulerHeun(); dt=dtmix, tstops=soloop.t, adaptive=false)

  @test - adj_soloop[end][1:length(u₀)] ≈ res_sde_u0  atol = 1e-14
  @test adj_soloop[end][length(u₀)+1:end] ≈ res_sde_p' atol = 1e-14

  adjprob = SDEAdjointProblem(sol,InterpolatingAdjoint(noisemixing=true, checkpointing=true),dg!,tarray, nothing)
  adj_sol = solve(adjprob,EulerHeun(); dt=dtmix, adaptive=false,tstops=soloop.t)

  @test adj_soloop[end] ≈ adj_sol[end]  rtol=1e-8

  adjprob = SDEAdjointProblem(sol,InterpolatingAdjoint(noisemixing=true, checkpointing=false),dg!,tarray, nothing)
  adj_sol = solve(adjprob,EulerHeun(); dt=dtmix, adaptive=false,tstops=soloop.t)

  @test adj_soloop[end] ≈  adj_sol[end]  rtol=1e-8
end

@testset "mutating non-diagonal noise" begin
  a!(du,u,_p,t) = (du .= -u)
  a(u,_p,t) = -u

  function b!(du,u,_p,t)
    KR, KI = _p[1:2]

    du[1,1] = KR
    du[2,1] = KI
  end

  function b(u,_p,t)
    KR, KI = _p[1:2]

    [ KR zero(KR)
      KI zero(KR) ]
  end

  p = [1.,0.]

  prob! = SDEProblem{true}(a!,b!,[0.,0.],(0.0,0.1),p,noise_rate_prototype=eltype(p).(zeros(2,2)))
  prob = SDEProblem{false}(a,b,[0.,0.],(0.0,0.1),p,noise_rate_prototype=eltype(p).(zeros(2,2)))

  function loss(p;SDEprob=prob,sensealg=BacksolveAdjoint())
    _prob = remake(SDEprob,p=p)
    sol = solve(_prob, EulerHeun(), dt=1e-5, sensealg=sensealg)
    return sum(Array(sol))
  end

  function compute_dp(p, SDEprob, sensealg)
    Random.seed!(seed)
    Zygote.gradient(p->loss(p,SDEprob=SDEprob,sensealg=sensealg), p)[1]
  end

  # test mutating against non-mutating

  # non-mutating

  dp1 = compute_dp(p, prob, ForwardDiffSensitivity())
  dp2 = compute_dp(p, prob, BacksolveAdjoint())
  dp3 = compute_dp(p, prob, InterpolatingAdjoint())

  @show dp1 dp2 dp3

  # different vjp choice
  _dp2 = compute_dp(p, prob, BacksolveAdjoint(autojacvec=ReverseDiffVJP()))
  @test dp2 ≈ _dp2 rtol=1e-8
  _dp3 = compute_dp(p, prob, InterpolatingAdjoint(autojacvec=ReverseDiffVJP()))
  @test dp3 ≈ _dp3 rtol=1e-8

  # mutating
  _dp1 = compute_dp(p, prob!, ForwardDiffSensitivity())
  _dp2 = compute_dp(p, prob!, BacksolveAdjoint())
  _dp3 = compute_dp(p, prob!, InterpolatingAdjoint())

  @test dp1 ≈ _dp1 rtol=1e-8
  @test dp2 ≈ _dp2 rtol=1e-8
  @test dp3 ≈ _dp3 rtol=1e-8
end


text:
    export isless

    """
    `Basics.interpolateOnGridGrasp92((F::Array{Float64,1}, oldGrid::Radial.Grid), (G::Void, newGrid::Radial.Grid) )`
        ... to interpolate the (radial) function F from the oldgrid to newgrid, by a call to the Grasp92 Fortran procedures;
            an function G::Array{Float64,1} is returned. **Not yet implemented !**
    """
    function Basics.interpolateOnGridGrasp92(from::Tuple{Array{Float64,1},Radial.Grid}, to::Tuple{DataType,Radial.Grid})
        # First prepare all fields and settings to call quad_grasp92()
        error("Not yet implemented !")
    end

    """
    `Basics.interpolateOnGridTrapezRule((F::Array{Float64,1}, oldGrid::Radial.Grid), (G::Void, newGrid::Radial.Grid)`
        ... to integrate the (radial) function F from the oldgrid to newgrid, by using a trapez rule from the Numerical Recipies.
            A function G::Array{Float64,1} is returned. **Not yet implemented !**.
    """
    function Basics.interpolateOnGridTrapezRule(from::Tuple{Array{Float64,1},Radial.Grid}, to::Tuple{DataType,Radial.Grid})
        error("Not yet implemented !")
    end

    """
    `Basics.integrate()`
        ... integrates a one- or two-dimensional function by different numerical methods, either on a given grid or for a general function.

    + `("function: on radial grid, Newton-Cotes", F::Array{Float64,1}, grid::Radial.Grid)`
        ... to integrate the (radial) function F over the given (radial) grid, int_0^infinity dr F(r) by using a 5-point Newton-Cotes
            integration formula; a value::Float64 is returned.

    + `("function: on radial grid, Simpson rule", F::Array{Float64,1}, grid::Radial.Grid)`
        ... to integrate the (radial) function F over the given (radial) grid, int_0^infinity dr F(r) by using Simpson's rule;
            a value::Float64 is returned.

    + `("function: on radial grid, trapez rule", F::Array{Float64,1}, grid::Radial.Grid)`
        ... to integrate the (radial) function F over the given (radial) grid, int_0^infinity dr F(r), by using a simple trapez rule;
            a value::Float64 is returned.
    """
    function Basics.integrate(sa::String, F::Array{Float64,1}, grid::Radial.Grid)
        if       sa == "function: on radial grid, Newton-Cotes"     wa = integrateOnGridNewtonCotes(F, grid)
        elseif   sa == "function: on radial grid, Simpson rule"     wa = integrateOnGridSimpsonRule(F, grid)
        elseif   sa == "function: on radial grid, trapez rule"      wa = integrateOnGridTrapezRule( F, grid)
        else     error("Unsupported keystring = $sa.")
        end

    end

    """
    `Basics.integrateOnGridNewtonCotes()`  ... integrates by using a 5-point Newton-Cotes integration formula; see Basic.integrate().
    """
    function Basics.integrateOnGridNewtonCotes(F::Array{Float64,1}, grid::Radial.Grid)
        coefficients = Array{Float64}([2 * 7, 32, 12, 32, 7]) * 2 / 45 * grid.h
        n = size(coefficients, 1)

        result = 0.;    i0     = 1

        while  abs(F[i0]) == 0
            i0 += 1
        end
        gamma  = log(F[i0]/F[i0+1] * grid.rp[i0+1]/grid.rp[i0]) / log(grid.r[i0]/grid.r[i0+1])
        result = 1/(gamma + 1) * grid.r[i0] * F[i0] / grid.rp[i0]

        result += F[i0] * coefficients[n]
        for i = i0+1:size(F, 1)
            result += coefficients[(i - i0) % (n - 1) + 1] * F[i]
        end

        if (size(F, 1) - i0) % (n - 1) + 1 == 1
            result -= coefficients[n] * F[size(F, 1)]
        end

        return( result )
    end

    """
    `Basics.integrateOnGridSimpsonRule()`  ... integrates by using Simpson's rule; see Basics.integrate().
    """
    function Basics.integrateOnGridSimpsonRule(F::Array{Float64,1}, grid::Radial.Grid)
        coefficients = Array{Float64}([2 * 1, 4, 1]) / 3 * grid.h
        n = size(coefficients, 1)

        result = 0.
        result += F[1] * coefficients[n]
        for i = 2:size(F, 1)
            result += coefficients[(i - 1) % (n - 1) + 1] * F[i]
        end

        if (size(F, 1) - 1) % (n - 1) + 1 == 1
            result -= coefficients[n] * F[size(F, 1)]
        end

        return( result )
    end

    """
    `Basics.integrateOnGridTrapezRule()`  ... integrates by using a simple trapez rule; see Basics.integrate().
    """
    function Basics.integrateOnGridTrapezRule(F::Array{Float64,1}, grid::Radial.Grid)
        coefficients = Array{Float64}([2 * 1, 1]) / 2 * grid.h
        n = size(coefficients, 1)

        result = 0.

        result += F[1] * coefficients[n]
        for i = 2:size(F, 1)
            result += coefficients[(i - 1) % (n - 1) + 1] * F[i]
        end

        if (size(F, 1) - 1) % (n - 1) + 1 == 1
            result -= coefficients[n] * F[size(F, 1)]
        end

        return( result )
    end

    function Base.isless(x::ManyElectron.Level, y::ManyElectron.Level)
        return x.energy < y.energy
    end

    function Base.isless(x::Hfs.IJF_Level, y::Hfs.IJF_Level)
        return x.energy < y.energy
    end

    function Base.isless(x::Basics.AngularJ64, y::Basics.AngularJ64)
        return x.num < y.num
    end

    function Base.isless(x::Cascade.AbsorptionCrossSection, y::Cascade.AbsorptionCrossSection)
        return x.photonEnergy < y.photonEnergy
    end

    """
    `Basics.isSimilar()`  ... returns true if two instances are similar to each other, and false otherwise.

    + `(keya::LevelKey, keyb::LevelKey, relAcc::Float64)`
        ... returns true if two level keys refer to the same level, i.e. level with the same symmetry and
            if the relative energy abs( (E_a - E_b)/E_a ) < relAcc. It returns false otherwise.
    """
    function Basics.isSimilar(keya::LevelKey, keyb::LevelKey, relAcc::Float64)
        if  keya.sym == keyb.sym   &&   abs( (keya.energy - keyb.energy)/keya.energy ) < relAcc    return(true)
        else                                                                                return(false)
        end
    end

    """
    `Basics.isSymmetric(matrix::Array{Float64,2})`
        ... returns true if the matrix is symmetric and false otherwise; if not symmetric, the function prints the
            first 10 pairs  M_ij, M_ji  which violate the symmetry condition.
    """
    function Basics.isSymmetric(matrix::Array{Float64,2})
        symmetric = true;    na = 0
        if  size(matrix,1) !=  size(matrix,2)  error("not quadratic")   end
        #
        for j = 1:size(matrix,2)
            for i = j+1:size(matrix,1)
                if abs( matrix[i,j] - matrix[j,i]) > 1.0e-10
                    na = na + 1;   symmetric = false;   println("*** $i  $j  $(matrix[i,j])  $(matrix[j,i]) ")
                    if  na >= 10   return( symmetric )  end
                end
            end
        end

        return( symmetric )
    end

    """
    `Basics.isZero(matrix::Array{Float64,2})`
        ... returns true if the matrix is zero and false otherwise; it sums up the |M[i,j]^2 and prints the results.
            If not zero, the function prints the first 10 matrix elements  M_ij  which violate this condition.
    """
    function Basics.isZero(matrix::Array{Float64,2})
        isz = true;    na = 0;   wa = 0.
        #
        for j = 1:size(matrix,2)
            for i = 1:size(matrix,1)
                wa = wa + abs( matrix[i,j])^2
                #
                if abs( matrix[i,j]) > 1.0e-10
                    na = na + 1;   isz = false;   println("*** $i  $j  $(matrix[i,j]) ")
                    if  na >= 10   return( isz )    end
                end
            end
        end
        println(">> isZero():   sum^2 = $(wa)")

        return( isz )
    end

    """
    `Basics.isViolated()`  ... returns true if some rule/limitation is violated, and false otherwise.

    + `(conf::Configuration, restriction::AbstractConfigurationRestriction)`
        ... returns true if restriction is violated by the given conf, and false otherwise.
    """
    function Basics.isViolated(conf::Configuration, restriction::AbstractConfigurationRestriction)
        wa = false;  ne = 0
        if      typeof(restriction) == RestrictNoElectronsTo
            for (sh,v) in  conf.shells
                if  sh.n >= restriction.nmin    ||    sh.l >= restriction.lmin    ne = ne + v   end
            end
            if  ne > restriction.ne                                   wa = true   end
        elseif  typeof(restriction) == RestrictParity
            if  Basics.determineParity(conf) != restriction.parity    wa = true   end
        elseif  typeof(restriction) == RestrictToShellDoubles
            for (sh,v) in  conf.shells
                if  sh.n >= restriction.nmin
                    if  sh.l >= restriction.lmin
                        println("restrictions:  $sh  $(sh.n)  $v  $(restriction.nmin)")
                        if !(v in [0,2])                              wa = true   end
                    end
                end
            end
        elseif  typeof(restriction) == RequestMinimumOccupation
            for (sh,v) in  conf.shells
                if  sh in restriction.shells   ne = ne + v     end
            end
            if  ne < restriction.ne                                   wa = true   end
        elseif  typeof(restriction) == RequestMaximumOccupation
            for (sh,v) in  conf.shells
                if  sh in restriction.shells   ne = ne + v     end
            end
            if  ne > restriction.ne                                   wa = true   end
        else    error("stop a")
        end

        return( wa )
    end

    """
    `Basics.isStandardSubshellList(basis::Basis)`
        ... returns true if the subshell list basis.subshells is in standard order, and false otherwise.
            It requests especially that both subshells of the same shell (n,l) occur in the sequence j = l-1/2, j = l+1/2
            and that l increases before n increases.
    """
    function Basics.isStandardSubshellList(basis::Basis)
        function  decimal(n,l,jnum)
            return( 10000n + 100l + jnum)
        end
        subshells = basis.subshells
        na = subshells[1].n;    la = Basics.subshell_l(subshells[1]);    ja = Basics.subshell_j(subshells[1])

        for  i = 2:length(subshells)
            nb = subshells[i].n;    lb = Basics.subshell_l(subshells[i]);    jb = Basics.subshell_j(subshells[1])
            if  decimal(nb,lb,jb.num) <= decimal(na,la,ja.num)    return( false )   end
        end

        return( true )
    end

    """
    Basics.merge(bases::Array{Basis,1})`
        ... to merge two (or more) atomic bases into a single basis::Basis. This method assumes the same number of electrons in all basis and
            that the subshell lists are the same or can be made 'consistent' to each other. Two bases have a consistent subshell list if all
            subshells, that appear in any of the two lists appear always in the same sequence (if they are not missing at all). In the merged
            basis, the radial orbitals are taken from the basis (in the bases-array}, from where they are found first.
    """
    function Basics.merge(bases::Array{Basis,1})

        if  length(bases) > 2
            println("Number of bases = $(length(bases)) ...")
            bs       = Basics.merge([ bases[1], bases[2] ])
            basesNew = [bs]
            for i = 3:length(bases)    push!( basesNew, bases[i])    end
            Basics.merge(basesNew)
        elseif  length(bases) == 1    return( basis[1] )
        elseif  length(bases) == 2
            # This is the essential step here; first create a consistent subshell list
            basisA = bases[1];    basisB = bases[2]
            wa = deepcopy( basisA.subshells );   wb = [ 100i  for i = 1:length(wa) ]
            # Add the subshells of basisB
            for  i = 1:length( basisB.subshells )
                if    basisB.subshells[i] in wa
                else  push!( wa, basisB.subshells[i] );    push!( wb, wb[i-1]+1 )
                end
            end
            # Now sort wb and create a new subshell list from wa due to this sorting
            wc = sortperm( wb )
            newSubshells = Subshell[]
            for  i in wc   push!( newSubshells, wa[i] )    end
            println("Sorted subshell list: $(newSubshells) ")
            error("Not yet fully implemented, see source code ")
            # Transform all CSF from basisA and basisB to newSubshells ... and apply CsfRExcludeDouble()
            # return new basis
            #
        else    error("stop a")
        end
    end

    """
    `Basics.merge(multiplets::Array{Multiplet,1})`
        ... to merge two (or more) atomic multiplets into a single multiplet::Multiplet. This method assumes (and checks) that all
            levels have level.hasStateRep = true and that all levels refer to the same basis.
    """
    function Basics.merge(multiplets::Array{Multiplet,1})

    if      length(multiplets) == 1    return( multiplets[1] )
        elseif  length(multiplets) == 2
            levels = copy(multiplets[1].levels);   basis = multiplets[1].levels[1].basis
            for  lev  in  multiplets[2].levels
                if  !(lev.hasStateRep)   ||   lev.basis != basis    error("Levels of multiplets do not refer to the same basis.")   end
                push!(levels, lev)
            end
            mp = Multiplet( multiplets[1].name * "+" * multiplets[2].name, levels )
            return( mp )
        elseif  length(multiplets) > 2
            mp = Basics.merge([ multiplets[1], multiplets[2] ])
            for  k = 3:length(multiplets)
                mp = Basics.merge([ mp, multiplets[k] ])
            end
            return( mp )
        else    error("stop a")
        end
    end

    """
    `Basics.modifyLevelEnergies(multiplet::Multiplet)`
        ... to shift the total energies of the atomic levels in multiplet interactively; a (new) multiplet::Multiplet is
            returned in which the total energies are just modified accordingly. **Not yet implemented !**
    """
    function Basics.modifyLevelEnergies(multiplet::Multiplet)
        error("Not yet implemented !")

        return( nothing )
    end

    # using Plots
    # pyplot()

    using RecipesBase

    """
    `Basics.plot()`  ... plots various quantities, often in a new window.

    + `("radial potentials", potentials::Array{Radial.Potential,1}, grid::Radial.Grid; N::Int64 = 0)`
        ... to plot one or more radial potentials, and where N::Int64 describes the number of grid points to be considered.
            call:  using Plots; pyplot()    ... to access this method by plot(...)
    """
    function Basics.plot(sa::String, potentials::Array{Radial.Potential,1}, grid::Radial.Grid; N::Int64 = 0)
        error("call instead:  using Plots; pyplot()    ... to access this method simply by plot(...)")
    end

    @recipe function f(sa::String, potentials::Array{Radial.Potential,1}, grid::Radial.Grid; N = 0)
        !(sa == "radial potentials")   &&   error("Unsupported keystring = $sa")
        wa = Float64[];   wc = [NaN for i=1:N];   labels = String[];   np = length(potentials)
        for  pot in potentials
            wb = wc
            nx = min(length(pot.Zr), N);   wb[1:nx] = pot.Zr[1:nx]
            append!(wa, wb)
            push!(labels, pot.name)
        end
        x = grid.r[1:N];     y = reshape(wa, (N, np))
        label --> reshape(labels, (1,np))
        x, y
    end

    """
    + `("radial orbitals: large", orbitals::Array{Radial.Orbital,1}, grid::Radial.Grid; N::Int64 = 0)`
        ... to plot the large component of one or more radial orbitals, and where N::Int64 describes the number of grid
            points to be considered.
    + `("radial orbitals: small", orbitals::Array{Radial.Orbital,1}, grid::Radial.Grid; N::Int64 = 0)`
        ... to plot the small component of one or more radial orbitals.
    + `("radial orbitals: both",  orbitals::Array{Radial.Orbital,1}, grid::Radial.Grid; N::Int64 = 0)`
        ... to plot the large and small component of one or more radial orbitals.

        call:  using Plots; pyplot()    ... to access this method by plot(...)
    """
    function Basics.plot(sa::String, orbitals::Array{Radial.Orbital,1}, grid::Radial.Grid; N = 0)
        error("call instead:  using Plots; pyplot()    ... to access this method simply by plot(...)")
    end

    @recipe function f(sa::String, orbitals::Array{Radial.Orbital,1}, grid::Radial.Grid; N = 0)
        wa = Float64[];   wc = [NaN for i=1:N];   labels = String[];   np = length(orbitals)

        if       sa == "radial orbitals: large"
            for  orb in orbitals
                wb = wc
                nx = min(length(orb.P), N);   wb[1:nx] = orb.P[1:nx]
                append!(wa, wb)
                push!(labels, "$(orb.subshell):large")
            end
            x = grid.r[1:N];     y = reshape(wa, (N, np))

        elseif   sa == "radial orbitals: small"
            for  orb in orbitals
                wb = wc
                nx = min(length(orb.Q), N);   wb[1:nx] = orb.Q[1:nx]
                append!(wa, wb)
                push!(labels, "$(orb.subshell):small")
            end
            x = grid.r[1:N];     y = reshape(wa, (N, np))

        elseif   sa == "radial orbitals: both"
            for  orb in orbitals
                wb = wc;    nx = min(length(orb.P), N);   wb[1:nx] = orb.P[1:nx];    append!(wa, wb)
                wb = wc;    nx = min(length(orb.Q), N);   wb[1:nx] = orb.Q[1:nx];    append!(wa, wb)
                push!(labels, "$(orb.subshell):large");   push!(labels, "$(orb.subshell):small")
            end
             x = grid.r[1:N];     y = reshape(wa, (N, 2np))

        else   error("Unsupported keystring = $sa")
        end
        label --> permutedims(labels)
        x, y
    end

    """
    + `("spectrum: transition rates over energy", lines::Array{PhotoEmission.Line,1})`
        ... to plot the transition rates of all lines as function of their transition energies. The plot is shown in a new
            window but nothing is returned otherwise. **Not yet implemented !**

    + `("spectrum: oscillator strength over energy, emission", lines::Array{PhotoEmission.Line,1})` or
        `("spectrum: oscillator strength over energy, absorption", lines::Array{PhotoEmission.Line,1})`
        ... to plot the absorption oscillator strength of all lines as function of their transition energies. Again, a new
            window is opened but nothing returned by this method. **Not yet implemented !**
    """
    function Basics.plot(sa::String, lines::Array{PhotoEmission.Line,1})
        error("call instead:  using Plots; pyplot()    ... to access this method simply by plot(...) ... not yet implemented !")
    end

    """
    + `("spectrum: transition rates over energy, Gaussian", lines::Array{PhotoEmission.Line,1}; widths=value::Float64)` or
        `("spectrum: transition rates over energy, Lorentzian", lines::Array{PhotoEmission.Line,1}; widths=value::Float64)`
        ... to plot the transition rates of all lines as function of their transition energies but with a Gaussian or Lorentzian
            distribution. Again, a new window is opened but nothing returned by this method. It still need to be decided how
            the widths (and, perhaps, other parameters) will be communicated to the method. **Not yet implemented !**
    """
    function Basics.plot(sa::String, lines::Array{PhotoEmission.Line,1}, widths::Float64)
        error("call instead:  using Plots; pyplot()    ... to access this method simply by plot(...) ... not yet implemented !")
    end

ext: jl


text:
# A simple spiking network simulator.

module SpikeNet

include("elemwise.jl")

abstract type Pathway end
export Pathway

function record! end
record!(r::Void, step) = nothing
export record!

function reset! end
export reset!

include("population.jl")
export Group

include("densepathway.jl")
export DensePathway

include("sparsepathway.jl")
export SparsePathway

include("spikerecorder.jl")
export RecordedSpikes

include("staterecorder.jl")
export RecordedState

include("correcorder.jl")
export JointStatsRecorder, covariance, variances, corrcoef, marginal_averages, joint_average

include("plotting.jl")
export raster_spikes

end


text:
using Revise
using Distributions, Random
using Flux, Flux.Tracker
using RCall

@rimport graphics as rgraphics
@rimport rcommon

compute_ke(x) = x^2/2

function hmc(curr_state::Real, log_prob::Function, L::Integer, eps::Real)
  U(q) = -log_prob(q)
  grad_U(q) = Tracker.gradient(U, q)[1]

  q_history = zeros(L + 1)
  p_history = zeros(L + 1)

  q = curr_state
  p = randn()

  q_history[1] = q
  p_history[1] = p

  curr_K = compute_ke(p)
  curr_U = U(q)

  p -= eps * Tracker.data(grad_U(q)) / 2
  for i in 1:L
    q += eps * p
    q_history[i + 1] = q
    if i < L
      p -= eps * Tracker.data(grad_U(q))
      p_history[i + 1] = p
    end
  end
  p -= eps * Tracker.data(grad_U(q)) / 2
  p_history[L + 1] = p

  cand_U = U(q)
  cand_K = compute_ke(p)

  log_accept_ratio = curr_U + curr_K - cand_U - cand_K

  if log_accept_ratio > log(rand())
    q, lp = q, -cand_U
  else
    q, lp = curr_state, -curr_U
  end

  return Dict(:q => q, :log_prob => lp,
              :qhist => q_history, :phist => p_history)
end

### MAIN ###

function simulate(; L, eps, n, sd=.5, init=0.0, plot_bound=4)
  log_prob(x) = logpdf(Normal(0, sd), x)
  rgraphics.plot(0, xlab="q", ylab="p", typ="n",
                 ylim=[-1, 1]*plot_bound, xlim=[-sd, sd]*plot_bound)
  sim = hmc(init, log_prob, L, eps)
  qs = zeros(n + 1)
  qs[1] = sim[:q]

  for i in 1:n
    sim = hmc(sim[:q], log_prob, L, eps)
    rgraphics.lines(sim[:qhist], sim[:phist], typ="l")
    rgraphics.points(sim[:qhist][end], sim[:phist][end], pch=4, lwd=2)
    qs[i + 1] = sim[:q]
  end

  return qs
end

qs = simulate(L=2^4, eps=.1, n=30, sd=.5)
rcommon.plotPost(qs, typ="l", xlab="", ylab="", main="");

function simulate_gamma(; L, eps, n, sd=.5, init=0.0, plot_bound=4)
  shape = 2
  dist = Gamma(shape, sd/sqrt(shape))
  log_prob(logx) = logpdf(dist, exp(logx)) + logx
  rgraphics.par(mfrow=[2, 1])
  x = collect(range(0, sd*plot_bound, step=.01))
  rgraphics.plot(x, pdf.(dist, x), ylab="", xlab="", main="", typ="l")
  rgraphics.plot(0, xlab="q", ylab="p", typ="n",
                 ylim=[-1, 1]*plot_bound, xlim=[0, sd]*plot_bound)
  sim = hmc(init, log_prob, L, eps)
  qs = zeros(n + 1)
  qs[1] = sim[:q]

  for i in 1:n
    sim = hmc(sim[:q], log_prob, L, eps)
    rgraphics.lines(exp.(sim[:qhist]), sim[:phist], typ="l");
    rgraphics.points(exp(sim[:qhist][end]), sim[:phist][end], pch=4, lwd=2);
    # rgraphics.abline(v=[exp(sim[:qhist][end]), sim[:phist][end]], lty=2, lwd=.5)
    # rgraphics.points(exp(sim[:qhist][1]), sim[:phist][1], pch="$i", lwd=2)
    qs[i + 1] = sim[:q]
  end
  rgraphics.abline(v=mean(dist), lty=2, col="grey")
  rgraphics.par(mfrow=[1, 1])

  return qs
end

Random.seed!(3);
qs = simulate_gamma(L=2^4, eps=.1, n=10, sd=1, plot_bound=5)
rcommon.plotPost(exp.(qs), typ="l", xlab="", ylab="", main="");
_ = nothing

# sim = hmc(0.0, logx -> logpdf(Gamma(2, 3), exp(logx)) + logx, 10, .1)
# rgraphics.plot(sim[:qhist], main="", xlab="", ylab="")


text:
module Evolutionary
    using Random, LinearAlgebra, Statistics
    using Base: @kwdef
    using UnPack: @unpack
    using NLSolversBase: NLSolversBase, AbstractObjective, ConstraintBounds,
                         AbstractConstraints, nconstraints_x, nconstraints

    import NLSolversBase: NonDifferentiable, f_calls, value, value!
    import Base: show, copy, minimum, summary, identity, getproperty, rand, getindex, length, copyto!, setindex!

    export AbstractStrategy, strategy, mutationwrapper,
           IsotropicStrategy, AnisotropicStrategy, NoStrategy,
           isfeasible, BoxConstraints, apply!, penalty, penalty!,
           PenaltyConstraints, WorstFitnessConstraints, MixedTypePenaltyConstraints,
           # ES mutations
           gaussian, cauchy,
           # GA mutations
           flip, bitinversion, domainrange, inversion, insertion, swap2, scramble, shifting, PM, MIPM,
           # ES recombinations
           average, marriage,
           # GA recombinations
           singlepoint, twopoint, uniform,
           discrete, waverage, intermediate, line, HX, LX, MILX,
           PMX, OX1, CX, OX2, POS,
           # GA selections
           ranklinear, uniformranking, roulette, rouletteinv, sus, susinv, tournament, truncation,
           # DE selections
           random, permutation, randomoffset, best,
           # DE recombinations
           uniformbin, exponential,
           # GP exports
           Terminal, subtree, point, hoist, shrink,
           # Optimization methods
           ES, CMAES, GA, DE, TreeGP, AFS
           # re-export
           NonDifferentiable, value, value!

    # optimize API
    include("api/types.jl")
    include("api/results.jl")
    include("api/utilities.jl")
    include("api/constraints.jl")
    include("api/optimize.jl")
    include("api/expressions.jl")

    # Evolution Strategy
    include("es.jl")
    include("cmaes.jl")

    # Genetic Algorithms
    include("ga.jl")

    # Differential Evolution
    include("de.jl")

    #Artificial Fish swarm
    include("afs.jl")

    # Genetic Programming
    include("api/protected.jl")
    include("gp.jl")

    # ES & GA recombination functions
    include("recombinations.jl")

    # ES & GA mutation functions
    include("mutations.jl")

    # GA selection functions
    include("selections.jl")

end


text:
using Jive
using Documenter

makedocs(
    build = joinpath(@__DIR__, "local" in ARGS ? "build_local" : "build"),
    modules = [Jive],
    clean = false,
    format = Documenter.HTML(prettyurls = !("local" in ARGS)),
    sitename = "Jive.jl 👣",
    authors = "WooKyoung Noh",
    pages = Any[
        "Home" => "index.md",
        "runtests" => "runtests.md",
        "watch" => "watch.md",
        "@skip" => "skip.md",
        "@onlyonce" => "onlyonce.md",
        "@If" => "If.md",
        "@useinside" => "useinside.md",
        "@mockup" => "mockup.md",
        "`@__END__`" => "END.md",
    ],
)


text:
export TrustRegion, Combi, decreaseFact, convert_TR, convert_TR!, extract

"Exception type raised in case of error."
mutable struct TrustRegionException <: Exception
    msg  :: String
end

mutable struct TrustRegion{T}
    α₀ :: T
    α :: T
    max_α :: T
    acceptance_threshold :: T
    increase_threshold :: T
    reduce_threshold :: T
    increase_factor :: T
    decrease_factor :: T
    max_unsuccinarow :: Int
    #params :: Tparams

    function TrustRegion(α₀ :: T;
                         max_α :: T = T(1.0)/ sqrt(eps(T)),
                         acceptance_threshold :: T = T(0.1),
                         increase_threshold :: T = T(0.75),
                         reduce_threshold :: T = T(0.1),
                         increase_factor :: T = T(5.0),
                         decrease_factor :: T = T(0.1),
                         max_unsuccinarow :: Int = 30) where T
                         #params :: Tparams = Void)

        α₀ > T(0.0) || (α₀ = T(1.0))
        max_α > α₀ || throw(TrustRegionException("Invalid α₀"))
        (T(0.0) < acceptance_threshold < increase_threshold < T(1.0)) || throw(TrustRegionException("Invalid thresholds"))
        (T(0.0) < decrease_factor < T(1.0) < increase_factor) || throw(TrustRegionException("Invalid decrease/increase factors"))

        return new{T}(α₀, α₀, max_α,
                   acceptance_threshold, increase_threshold, reduce_threshold,
                   increase_factor, decrease_factor, max_unsuccinarow)#, params)
    end
end

"""Compute the actual vs. predicted reduction radio ∆f/Δm, where
Δf = f - f_trial is the actual reduction is an objective/merit/penalty function,
Δq = q - q_trial is the reduction predicted by the model q of f.
We assume that q is being minimized, and therefore that Δq > 0.
"""
function compute_r(nlp, f, Δf, Δq, slope, d, xnext, gnext, robust)
    # If ever the next gradient is computed for round off errors reason, signal it.
    # # printstyled("dans compute r eltype(f) = $(eltype(f)) \n", color = :yellow)
    # # @show eltype(f)
    T = eltype(f)
    good_grad = false
    if robust & ((Δq < 10000 * eps(T)) | (abs(Δf) < 10000 * eps(T) * abs(f)) )
        # trap potential cancellation errors
        grad!(nlp, xnext, gnext);
        good_grad = true
        slope_next = dot(gnext, d)

        Δf = -(slope_next+slope) / 2.0
    end
    r = Δf / Δq
    if isnan(r)
        r = 0.0;
    end;

    return r, good_grad, gnext
end

# A few headers displays
function print_header(Title)

    @printf(" \n\n %-15s\n",Title)
    @printf("\n%-16s  %5s  %9s  %7s %7s  %5s  %5s  %5s  %6s  %s   %s\n",
            "Name", "nvar", "f", "‖∇f‖∞","‖∇f‖₂", "#obj", "#grad", " H  ","#hprod", "status", "time")
end

function display_header_problems()
    @printf("%-15s  %5s  %9s  %7s %7s  %5s  %5s  %6s  %s\n",
            "Name", "nvar", "f", "‖∇f‖∞","‖∇f‖₂", "#obj", "#grad", "#hprod", "status")
end

function display_header_iterations()
    @printf(" Iter    f          ||g||       λ        α       status \n",)
end

# A few displays for verbose iterations
function display_failure(iter, fnext, λ, α)
    @printf("%4d  %10.3e             %7.1e %7.1e    unsuccessful\n", iter, fnext, λ, α)
end

function display_v_success(iter, f, norm_g, λ, α)
    @printf("%4d  %10.3e %9.2e   %7.1e %7.1e Very successful\n", iter, f, norm_g, λ, α)
end

function display_success(iter, f, norm_g, λ, α)
    @printf("%4d  %10.3e %9.2e   %7.1e %7.1e      successful\n", iter, f, norm_g, λ, α)
end

function print_stats(prob, dim, f, gNorm, gnorm2, calls, status, timt)
    @printf("%-16s  %5d  %9.2e  %7.1e  %7.1e  %5d  %5d %5d %6d  %s  %8.3f\n",
            prob, dim, f, gNorm, gnorm2, calls[1], calls[2], calls[3], calls[4], status, timt)
end

# default increase and decrease functions.
function decreaseBase(α:: T, TR:: TrustRegion) where T
    return α * TR.decrease_factor
end

function increase(α:: T, TR:: TrustRegion) where T
    return min(α * TR.increase_factor, TR.max_α)
end

function decreaseGen(X :: TPData, α:: T, TR:: TrustRegion) where T
    return decreaseBase(α, TR)
end

function decreaseFact(X :: PDataFact, α:: T, TR:: TrustRegion) where T
    X.success = false
    return decreaseBase(α, TR)
end

function increase(X :: TPData, α:: T, TR:: TrustRegion) where T
    return increase(α, TR)
end

stop_norm(x) = norm(x, Inf)

# Valid combinations
#
mutable struct Combi{T}
    hessian_rep :: Function
    PData       :: DataType
    solve_model :: Function
    pre_process :: Function
    decrease    :: Function
    params      :: Union{Tparams{T}, Tparams}
end

function extract(c :: Combi)
    return c.hessian_rep, c.PData, c.solve_model, c.pre_process, c.decrease, c.params
end

function convert_TR(T, TR_init :: TrustRegion)
	max_α_T                = T(TR_init.max_α)
	acceptance_threshold_T = T(TR_init.acceptance_threshold)
	increase_threshold_T   = T(TR_init.increase_threshold)
	reduce_threshold_T     = T(TR_init.reduce_threshold)
	increase_threshold_T   = T(TR_init.increase_threshold)
	decrease_factor_T      = T(TR_init.decrease_factor)

	TRT = TrustRegion(T(TR_init.α), max_α = max_α_T, acceptance_threshold = acceptance_threshold_T, increase_threshold = increase_threshold_T, decrease_factor = decrease_factor_T)

	return TRT
end


text:
using Distributions, Random
Random.seed!(1234)

# Sample from conditional X|Y=y
# Inputs:
# - y: The current value of y (Y_k)
# - μ: The mean vector of the joint density
# - Σ: The covariance matrix of X,Y
# Outputs:
# X_{k+1} where X_{k+1} ∼ f(X|Y=y)
function condition_on_Y(y, μ, Σ)
    μX = μ[1]
    μY = μ[2]
    σX2, σXY, _, σY2 = Σ

    conditional_mean = μX + σXY/σY2*(y-μY)
    conditional_variance = (1-σXY^2/(σX2*σY2))*σX2

    distribution = Normal(conditional_mean, conditional_variance)

    return rand(distribution)
end

# Sample from conditional Y|X=y
# Inputs:
# - x: The current value of x (X_{k+1})
# - μ: The mean vector of the joint density
# - Σ: The covariance matrix of X,Y
# Outputs:
# Y_{k+1} where Y_{k+1} ∼ f(y|X=x)
function condition_on_X(x, μ, Σ)
    μX = μ[1]
    μY = μ[2]
    σX2, σXY, _, σY2 = Σ

    conditional_mean = μY + σXY/σX2*(x-μX)
    conditional_variance = (1-σXY^2/(σX2*σY2))*σY2

    distribution = Normal(conditional_mean, conditional_variance)

    return rand(distribution)
end

# Gibbs sampler to obtain samples from (X,Y) ~ N(μ,Σ)
# Inputs:
# - μ: The joint density mean vector
# - Σ: The covaraince matrix of X and Y
# - n: The number of samples to return
# Outputs:
# An (n x 2) matrix where each row represent a sample
function bivariate_gibbs(μ, Σ, n)
    # Store samples here
    S = zeros(n,2)

    # Initialize S_0 to means
    S[1,:] = μ

    for i=2:n
        X_k1 = condition_on_Y(S[i-1,2], μ, Σ) # Get X_{k+1}
        Y_k1 = condition_on_X(X_k1, μ, Σ) # Get Y_{k+1}
        S[i,:] = [X_k1, Y_k1]
    end

    return S
end

# Define joint density parameters
μX = -1
μY = 1
σX = √2
σY = √3
σXY = 2
μ = [μX, μY]
Σ = [σX^2 σXY; σXY σY^2]

# Get samples
samples = bivariate_gibbs(μ, Σ, 10000)

# Create 2d histogram
using Seaborn, Pandas
pygui(true)
df = DataFrame(Dict(:X=>samples[:,1], :Y=>samples[:,2]))
jointplot(x = "X", y = "Y", data = df, s = 1, alpha = 0.18)

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

@testset "mm" begin
    @testset for elty in [Float32,Float64,Complex64,Complex128]
        A = sparse(rand(elty,m,k))
        B = rand(elty,k,n)
        C = rand(elty,m,n)
        alpha = rand(elty)
        beta = rand(elty)
        @testset "csr" begin
            d_B = CudaArray(B)
            d_C = CudaArray(C)
            d_A = CudaSparseMatrixCSR(A)
            @test_throws DimensionMismatch CUSPARSE.mm!('T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_B,'O')
            d_D = CUSPARSE.mm('N',alpha,d_A,d_B,beta,d_C,'O')
            h_D = to_host(d_D)
            D = alpha * A * B + beta * C
            @test D ≈ h_D
            d_D = CUSPARSE.mm('N',d_A,d_B,beta,d_C,'O')
            h_D = to_host(d_D)
            D = A * B + beta * C
            @test D ≈ h_D
            d_D = CUSPARSE.mm('N',d_A,d_B,d_C,'O')
            h_D = to_host(d_D)
            D = A * B + C
            @test D ≈ h_D
            d_D = CUSPARSE.mm('N',alpha,d_A,d_B,'O')
            h_D = to_host(d_D)
            D = alpha * A * B
            @test D ≈ h_D
            d_D = CUSPARSE.mm('N',d_A,d_B,'O')
            h_D = to_host(d_D)
            D = A * B
            @test D ≈ h_D
            d_D = d_A*d_B
            h_D = to_host(d_D)
            D = A * B
            @test D ≈ h_D
            @test_throws DimensionMismatch CUSPARSE.mm('T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm('N',alpha,d_A,d_B,beta,d_B,'O')
        end
        @testset "csc" begin
            d_B = CudaArray(B)
            d_C = CudaArray(C)
            d_A = CudaSparseMatrixCSC(A)
            @test_throws DimensionMismatch CUSPARSE.mm!('T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_B,'O')
            d_D = CUSPARSE.mm('N',alpha,d_A,d_B,beta,d_C,'O')
            h_D = to_host(d_D)
            D = alpha * A * B + beta * C
            @test D ≈ h_D
            d_D = CUSPARSE.mm('N',d_A,d_B,beta,d_C,'O')
            h_D = to_host(d_D)
            D = A * B + beta * C
            @test D ≈ h_D
            d_D = CUSPARSE.mm('N',d_A,d_B,d_C,'O')
            h_D = to_host(d_D)
            D = A * B + C
            @test D ≈ h_D
            d_D = CUSPARSE.mm('N',alpha,d_A,d_B,'O')
            h_D = to_host(d_D)
            D = alpha * A * B
            @test D ≈ h_D
            d_D = CUSPARSE.mm('N',d_A,d_B,'O')
            h_D = to_host(d_D)
            D = A * B
            @test D ≈ h_D
            d_D = d_A*d_B
            h_D = to_host(d_D)
            D = A * B
            @test D ≈ h_D
            @test_throws DimensionMismatch CUSPARSE.mm('T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm('N',alpha,d_A,d_B,beta,d_B,'O')
        end
    end
end

@testset "mm_symm" begin
    @testset for elty in [Float32,Float64,Complex64,Complex128]
        A = sparse(rand(elty,m,m))
        A = A + A.'
        B = rand(elty,m,n)
        C = rand(elty,m,n)
        alpha = rand(elty)
        beta = rand(elty)
        @testset "csr" begin
            d_B = CudaArray(B)
            d_C = CudaArray(C)
            d_A = Symmetric(CudaSparseMatrixCSR(A))
            d_C = CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_C,'O')
            h_C = to_host(d_C)
            D = alpha * A * B + beta * C
            @test D ≈ h_C
            d_C = d_A.' * d_B
            h_C = to_host(d_C)
            D = A.' * B
            @test D ≈ h_C
            d_B = CudaArray(rand(elty,k,n))
            @test_throws DimensionMismatch CUSPARSE.mm!('T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_B,'O')
        end
        @testset "csc" begin
            d_B = CudaArray(B)
            d_C = CudaArray(C)
            d_A = Symmetric(CudaSparseMatrixCSC(A))
            d_C = CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_C,'O')
            h_C = to_host(d_C)
            D = alpha * A * B + beta * C
            @test D ≈ h_C
            d_C = d_A.' * d_B
            h_C = to_host(d_C)
            D = A.' * B
            @test D ≈ h_C
            d_B = CudaArray(rand(elty,k,n))
            @test_throws DimensionMismatch CUSPARSE.mm!('T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_B,'O')
        end
    end
end
@testset "mm_herm" begin
    @testset for elty in [Float32,Float64,Complex64,Complex128]
        A = sparse(rand(elty,m,m))
        A = A + A'
        B = rand(elty,m,n)
        C = rand(elty,m,n)
        alpha = rand(elty)
        beta = rand(elty)
        @testset "csr" begin
            d_B = CudaArray(B)
            d_C = CudaArray(C)
            d_A = Hermitian(CudaSparseMatrixCSR(A))
            d_C = CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_C,'O')
            h_C = to_host(d_C)
            D = alpha * A * B + beta * C
            @test D ≈ h_C
            d_B = CudaArray(rand(elty,k,n))
            @test_throws DimensionMismatch CUSPARSE.mm!('T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_B,'O')
        end
        @testset "csc" begin
            d_B = CudaArray(B)
            d_C = CudaArray(C)
            d_A = Hermitian(CudaSparseMatrixCSC(A))
            d_C = CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_C,'O')
            h_C = to_host(d_C)
            D = alpha * A * B + beta * C
            @test D ≈ h_C
            d_B = CudaArray(rand(elty,k,n))
            @test_throws DimensionMismatch CUSPARSE.mm!('T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_B,'O')
        end
    end
end

@testset "mm2" begin
    @testset for elty in [Float32,Float64,Complex64,Complex128]
        A = sparse(rand(elty,m,k))
        B = rand(elty,k,n)
        C = rand(elty,m,n)
        alpha = rand(elty)
        beta = rand(elty)
        @testset "csr" begin
            d_B = CudaArray(B)
            d_C = CudaArray(C)
            d_A = CudaSparseMatrixCSR(A)
            @test_throws DimensionMismatch CUSPARSE.mm2!('N','T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm2!('T','N',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm2!('T','T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm2!('N','N',alpha,d_A,d_B,beta,d_B,'O')
            d_D = CUSPARSE.mm2('N','N',alpha,d_A,d_B,beta,d_C,'O')
            h_D = to_host(d_D)
            D = alpha * A * B + beta * C
            @test D ≈ h_D
            d_D = CUSPARSE.mm2('N','N',d_A,d_B,beta,d_C,'O')
            h_D = to_host(d_D)
            D = A * B + beta * C
            @test D ≈ h_D
            d_D = CUSPARSE.mm2('N','N',d_A,d_B,d_C,'O')
            h_D = to_host(d_D)
            D = A * B + C
            @test D ≈ h_D
            d_D = CUSPARSE.mm2('N','N',alpha,d_A,d_B,'O')
            h_D = to_host(d_D)
            D = alpha * A * B
            @test D ≈ h_D
            d_D = CUSPARSE.mm2('N','N',d_A,d_B,'O')
            h_D = to_host(d_D)
            D = A * B
            @test D ≈ h_D
        end
        @testset "csc" begin
            d_B = CudaArray(B)
            d_C = CudaArray(C)
            d_A = CudaSparseMatrixCSC(A)
            @test_throws DimensionMismatch CUSPARSE.mm2!('N','T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm2!('T','N',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm2!('T','T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm2!('N','N',alpha,d_A,d_B,beta,d_B,'O')
            d_D = CUSPARSE.mm2('N','N',alpha,d_A,d_B,beta,d_C,'O')
            h_D = to_host(d_D)
            D = alpha * A * B + beta * C
            @test D ≈ h_D
            d_D = CUSPARSE.mm2('N','N',d_A,d_B,beta,d_C,'O')
            h_D = to_host(d_D)
            D = A * B + beta * C
            @test D ≈ h_D
            d_D = CUSPARSE.mm2('N','N',d_A,d_B,d_C,'O')
            h_D = to_host(d_D)
            D = A * B + C
            @test D ≈ h_D
            d_D = CUSPARSE.mm2('N','N',alpha,d_A,d_B,'O')
            h_D = to_host(d_D)
            D = alpha * A * B
            @test D ≈ h_D
            d_D = CUSPARSE.mm2('N','N',d_A,d_B,'O')
            h_D = to_host(d_D)
            D = A * B
            @test D ≈ h_D
        end
    end
end
@testset "bsrmm2" begin
    @testset for elty in [Float32,Float64,Complex64,Complex128]
        A = sparse(rand(elty,m,k))
        B = rand(elty,k,n)
        C = rand(elty,m,n)
        alpha = rand(elty)
        beta = rand(elty)
        d_B = CudaArray(B)
        d_C = CudaArray(C)
        d_A = CudaSparseMatrixCSR(A)
        d_A = CUSPARSE.switch2bsr(d_A,convert(Cint,blockdim))
        @test_throws DimensionMismatch CUSPARSE.mm2('N','T',alpha,d_A,d_B,beta,d_C,'O')
        @test_throws DimensionMismatch CUSPARSE.mm2('T','N',alpha,d_A,d_B,beta,d_C,'O')
        @test_throws DimensionMismatch CUSPARSE.mm2('T','T',alpha,d_A,d_B,beta,d_C,'O')
        @test_throws DimensionMismatch CUSPARSE.mm2('N','N',alpha,d_A,d_B,beta,d_B,'O')
        d_D = CUSPARSE.mm2('N','N',alpha,d_A,d_B,beta,d_C,'O')
        h_D = to_host(d_D)
        D = alpha * A * B + beta * C
        @test D ≈ h_D
        d_D = CUSPARSE.mm2('N','N',d_A,d_B,beta,d_C,'O')
        h_D = to_host(d_D)
        D = A * B + beta * C
        @test D ≈ h_D
        d_D = CUSPARSE.mm2('N','N',d_A,d_B,d_C,'O')
        h_D = to_host(d_D)
        D = A * B + C
        @test D ≈ h_D
        d_D = CUSPARSE.mm2('N','N',alpha,d_A,d_B,'O')
        h_D = to_host(d_D)
        D = alpha * A * B
        @test D ≈ h_D
        d_D = CUSPARSE.mm2('N','N',d_A,d_B,'O')
        h_D = to_host(d_D)
        D = A * B
        @test D ≈ h_D
    end
end


text:
"""

# `node_edges`

$(SIGNATURES)

Internal.
"""
function node_edges(p::Path, s::Symbol, l::Symbol; debug=false)
  edges = Vector{Symbol}[]
  for sym in names(p.ug)[1]
    if p.ug[s, sym] == 1 && !(sym == l)
      append!(edges, [[s, sym]])
    end
  end
  edges
end

"""

# `all_edges`

$(SIGNATURES)

Part of the API, exported.
"""
function all_paths(d::DAG, f::Symbol, l::Symbol; debug=false)
  paths = Vector{Symbol}[]
  p = Path(d, f, l)
  debug && println(p)
  stack = Path[]
  push!(stack, p)
  debug && println("Starting stack = $(stack)\n")
  while !isempty(stack)
    debug && println("stack = $(stack)\n")
    p = pop!(stack)
    debug && println("stack = $(stack)\n")
    if length(p.visited) == 1 && p.visited == p.next
      p.next = node_edges(p, p.next[1][2], p.next[1][1]; debug=debug)
      debug && println(length(p.next))
      if length(p.next) > 1
        # Add a new Path to the stack
        for next in p.next[2:end]
          debug && println(next)
          np = deepcopy(p)
          np.next = p.next
          push!(stack, np)
          debug && println(stack)
        end
      end
    end
    debug && println("p.next = $(p.next)\n")
    for edge in p.next
      setdiff!(p.next, [edge])
      tmp = node_edges(p, edge[2], edge[1])
      while !isempty(tmp)
        pnew = deepcopy(p)
        push!(pnew.visited, edge)       # Edge to visited
        push!(pnew.path, edge[2])       # Neighbor to path
        next_edge = pop!(tmp)
        setdiff!(pnew.next, [next_edge])
        if !(next_edge in pnew.visited)
          if next_edge[end] in pnew.path
            #println("Drop path, back to initial start.")
          elseif next_edge[end] == l
            push!(pnew.visited, [next_edge[2], l])
            push!(pnew.path, l)
            push!(paths, pnew.path)
          else
            pnew.next = [next_edge]
            push!(stack, pnew)
          end
        end
      end
    end
  end
  paths
end

export
  all_paths


text:
using Test
using Mods

@testset "Constructors" begin
    @test Mod{17}() == 0
    @test Mod{17}(1) == GaussMod{17}(1,0)
    @test Mod{17}(1,2) == 1 + 2im
    a = Mod{17}(3)
    @test typeof(a) == Mod{17}
    a = Mod{17}(3,2)
    @test typeof(a) == GaussMod{17}
    a = GaussMod{17}()
    @test typeof(a) == GaussMod{17}
end

@testset "Mod arithmetic" begin
    p = 23
    a = Mod{p}(2)
    b = Mod{p}(25)
    @test a == b
    @test a == 2
    @test a == -21

    b = Mod{p}(20)
    @test a + b == 22
    @test a - b == -18
    @test a + a == 2a
    @test 0 - a == -a

    @test a * b == Mod{p}(17)
    @test (a / b) * b == a
    @test (b // a) * (2 // 1) == b
    @test a * (2 // 3) == (2a) * inv(Mod{p}(3))

    @test is_invertible(a)
    @test !is_invertible(Mod{10}(4))

    @test a^(p - 1) == 1
    @test a^(-1) == inv(a)
end

@testset "GaussMod arithmetic" begin
    p = 23
    a = GaussMod{p}(3 - im)
    b = GaussMod{p}(5 + 5im)

    @test a + b == 8 + 4im
    @test a + Mod{p}(11) == Mod{p}(14, 22)
    @test -a == 20 + im
    @test a - b == Mod{p}(3 - im - 5 - 5im)

    @test a * b == Mod{p}((3 - im) * (5 + 5im))
    @test a / b == Mod{p}((3 - im) // (5 + 5im))

    @test a^(p * p - 1) == 1
    @test is_invertible(a)
    @test a * inv(a) == 1

    @test a / (1 + im) == a / GaussMod{p}(1 + im)
    @test imag(a * a') == 0

end

@testset "Large Modulus" begin

    p = 9223372036854775783   # This is a large prime
    x = Mod{p}(-2)
    @test x * x == 4
    @test x + x == -4
    @test x / x == 1
    @test x / 3 == x / Mod{p}(3)
    @test (x / 3) * (3 // x) == 1
    @test x // x == value(x) / x
    @test x^4 == 16
    @test x^(p - 1) == 1   # Fermat Little Theorem test
    @test 2x == x + x
    @test x - x == 0
    y = inv(x)
    @test x * y == 1
    @test x + p == x
    @test x * p == 0
    @test p - x == x - 2x

    p = 9223372036854775783   # This is a large prime
    x = Mod{p}(-2 + 0im)
    @test x * x == 4
    @test x + x == -4
    @test x / x == 1
    @test x / 3 == x / Mod{p}(3)
    @test (x / 3) * (3 // x) == 1
    @test x // x == value(x) / x
    @test x^4 == 16
    @test x^(p - 1) == 1   # Fermat Little Theorem test
    @test 2x == x + x
    @test x - x == 0
    y = inv(x)
    @test x * y == 1
    @test x + p == x
    @test x * p == 0
    @test p - x == x - 2x

    @test 0 <= value(rand(Mod{p})) < p

end

@testset "CRT" begin
    p = 23
    q = 91
    a = Mod{p}(17)
    b = Mod{q}(32)

    x = CRT(a, b)
    @test a == mod(value(x), p)
    @test b == mod(value(x), q)

    c = Mod{101}(86)
    x = CRT(a,b,c)

    @test a == mod(value(x), p)
    @test b == mod(value(x), q)
    @test c == mod(value(x), 101)

end

@testset "Matrices" begin
    M = zeros(Mod{11}, 3, 3)
    @test sum(M) == 0

    M = ones(Mod{11}, 5, 5)
    @test sum(M) == 3

    M = rand(GaussMod{11}, 5, 6)
    @test size(M) == (5, 6)

    A = rand(Mod{17}, 5, 5)
    X = values.(A)
    @test sum(X) == sum(Mod{17}.(A))
end

@testset "Hashing/Iterating" begin
    x = Mod{17}(11)
    y = x + 0im
    @test x == y
    @test hash(x) == hash(y)
    @test typeof(x) !== typeof(y)

    A = Set([x,y])
    @test length(A) == 1

    v = [Mod{10}(t) for t = 1:15]
    w = [Mod{10}(t + 0im) for t = 1:15]
    S = Set(v)
    T = Set(w)
    @test length(S) == 10
    @test S == T
    @test union(S, T) == intersect(S, T)
end


text:
@recipe function f(C::CorrData)
    # set up the subplots
    grid := false
    layout := (2,1)

    lags = -C.maxlag:1/C.fs:C.maxlag
    times = Dates.format.(Dates.unix2datetime.(C.t),"yyyy/m/d HH:MM")

    # main heatmap
    @series begin
        legend := false
        seriestype := :heatmap
        ytickfontsize --> 10
        xtickfontsize --> 8
        subplot := 1
        seriescolor --> :balance
        lags,times,C.corr'
    end

    # bottom stack
    @series begin
        subplot := 2
        label := false
        title --> C.name
        xguide --> "Lag [s]"
        titlefontsize --> 6
        ytickfontsize --> 10
        xtickfontsize --> 8
        seriescolor --> :black
        xlims --> (lags[1],lags[end])
        lags,stack(C,allstack=true).corr
    end
end


text:
using Test
using DirectedAcyclicGraphs: map_values, groupby

@testset "Helper Functions" begin

    d = Dict("a" => 1, "b" => 3)
    @test map_values(x -> x+1, d, Int)["b"] == 4

    @test groupby(isodd, [1,2,3,4,5])[true] == [1,3,5]

end


text:
export MFABC, Particle, Cloud, BenchmarkCloud, MFABCParticle, MFABCCloud
export get_eta, phi, MakeMFABCCloud, MakeBenchmarkCloud, cost, length
export write_cloud

##### Everything below assumes exactly two fidelities:
# Future work will adapt methodology to a true multifidelity approach.

struct MFABC
    parameter_sampler::Function # Parameter sampler
    lofi::Function              # Map parameter to distance from observed data (low fidelity model) and coupling output
    hifi::Function              # Map parameter and coupling input to distance from observed data (high fidelity model)
end

struct Particle{N}
    k::Parameters
    dist::NTuple{N,Float64}
    cost::NTuple{N,Float64}
end
struct MFABCParticle
    p::Particle
    eta::Float64
    w::Float64
end

BenchmarkCloud = Array{Particle{2}, 1}
MFABCCloud = Array{MFABCParticle, 1}
Cloud = Union{BenchmarkCloud, MFABCCloud}

import Base.length
length(p::Particle) = length(p.cost)
length(pp::MFABCParticle) = length(pp.p)
length(c::Cloud, i::Integer) = length(filter(p->(length(p)==i), c))
cost(p::Particle) = sum(p.cost)
cost(p::Particle, i::Integer) = (i<=length(p) ? p.cost[i] : 0)
cost(pp::MFABCParticle) = cost(pp.p)
cost(pp::MFABCParticle, i) = cost(pp.p, i)
cost(c::Cloud) = sum(cost.(c))
cost(c::Cloud, i) = sum(cost.(c, i))

accept(p::Particle, epsilon::Float64) = (p.dist[2]<=epsilon)
accept_rate(c::BenchmarkCloud, epsilon::Float64) = count(p->accept(p,epsilon), c)/length(c)

######## Running simulations: benchmark particles ignore MFABC
function BenchmarkParticle(mfabc::MFABC, i::Int64=1)::Particle{2}

    k::Parameters = mfabc.parameter_sampler()
    (d_lo,pass),c_lo = @timed mfabc.lofi(k)
    (d_hi),c_hi = @timed mfabc.hifi(k,pass)

    return Particle{2}(k, (d_lo, d_hi), (c_lo, c_hi))
end

using Distributed
function MakeBenchmarkCloud(mfabc::MFABC, N::Int64)::BenchmarkCloud
# Simulate a Benchmark cloud for fixed continuation probabilities
    if nworkers()>1
        return pmap(i->BenchmarkParticle(mfabc,i), 1:N)
    else
        return map(i->BenchmarkParticle(mfabc,i), 1:N)
    end
end
function MakeBenchmarkCloud(mfabc::MFABC, N::Int64, fn::String)::BenchmarkCloud
    output = MakeBenchmarkCloud(mfabc, N)
    write_cloud(output, fn)
    return output
end

######## Running simulations: apply MFABC with known eta
# Apply either to:
# de novo simulations (mfabc::MFABC) or;
# previously completed simulations (p::Particle{2})

############
# PARTICLES
############

# De novo
function MFABCParticle(mfabc::MFABC, epsilons::Tuple{Float64, Float64}, etas::Tuple{Float64, Float64}, i::Int64=1)
    k::Parameters = mfabc.parameter_sampler()
    (d_lo,pass),c_lo = @timed mfabc.lofi(k)
    eta, w = (d_lo < epsilons[1]) ? (etas[1], 1) : (etas[2], 0)

    if rand()<eta
        (d_hi), c_hi = @timed mfabc.hifi(k,pass)
        close = ((d_lo, d_hi) .< epsilons)
        if xor(close...)
            w += (close[2]-close[1])/eta
        end
        p = Particle(k, (d_lo,d_hi), (c_lo,c_hi))
    else
        p = Particle(k, (d_lo,), (c_lo,))
    end
    return MFABCParticle(p,eta,w), sum(p.cost)
end

# Convert
function MFABCParticle(p::Particle{2}, epsilons::Tuple{Float64, Float64}, etas::Tuple{Float64, Float64}, i::Int64=1)
    close = (p.dist.<epsilons)
    etas==(1.0, 1.0) && (return MFABCParticle(p, 1.0, 1.0*close[2]), sum(p.cost))
    eta, w = (close[1]) ? (etas[1], 1.0) : (etas[2], 0.0)
    if rand()<eta
        if xor(close...)
            w += (close[2]-close[1])/eta
        end
        return MFABCParticle(p, eta, w), sum(p.cost)
    else
        q = Particle{1}(p.k, (p.dist[1],), (p.cost[1],))
        return MFABCParticle(q, eta, w), p.cost[1]
    end
end

############
# CLOUDS
############

## Known, fixed eta (no dependence on a function to be estimated, either)

# de novo
function MakeMFABCCloud(mfabc::MFABC, epsilons::Tuple{Float64, Float64}, etas::Tuple{Float64, Float64}; kwargs...)::MFABCCloud
# Multifidelity ABC algorithm with known continuation probabilities (which can therefore be parallelised)
    kwargs = Dict(kwargs)
    if haskey(kwargs, :N)
        N = Integer(kwargs[:N]) # Specifies the number of particles to simulate
        if nworkers()>1
            return pmap(i->MFABCParticle(mfabc, epsilons, etas, i)[1], 1:N)
        else
            return map(i->MFABCParticle(mfabc, epsilons, etas, i)[1], 1:N)
        end
    elseif haskey(kwargs, :budget)
        b = Float64(kwargs[:budget])
        run_cost = 0.0
        cloud = MFABCCloud()
        while 0<1
            pp, c = MFABCParticle(mfabc, epsilons, etas)
            run_cost += c
            (run_cost < b) ? push!(cloud, pp) : return cloud
        end
    else
        error("No size indication (N or budget)")
    end
end

# convert
function MakeMFABCCloud(s::BenchmarkCloud, epsilons::Tuple{Float64, Float64}, etas::Tuple{Float64, Float64})::MFABCCloud
    return map(p->MFABCParticle(p, epsilons, etas)[1], s) # First component is the particle (second is cost of particle)
end
function MakeMFABCCloud(s::BenchmarkCloud, epsilons::Tuple{Float64, Float64}, etas::Tuple{Float64, Float64}, budget::Float64)::MFABCCloud
    function truncate(cld, b, args...)
        return cld[1:searchsortedlast(cumsum(cost.(cld, args...)), b)]
    end
    mf = MakeMFABCCloud(truncate(s,budget,1), epsilons, etas)
    return truncate(mf, budget)
end

## Unknown eta: need to take into account finding a good eta, potentially given a function to be estimated

# de novo: this has to create a burn-in cloud and then adapt eta afterwards
function MakeMFABCCloud(mfabc::MFABC, epsilons::Tuple{Float64, Float64}; burnin::Int64, method::String="mf", kwargs...)::MFABCCloud
    s = MakeBenchmarkCloud(mfabc, burnin)
    return MakeMFABCCloud(mfabc, s, epsilons; method=method, kwargs...)
end
function MakeMFABCCloud(mfabc::MFABC, s::BenchmarkCloud, epsilons::Tuple{Float64, Float64}; method::String="mf", kwargs...)
# Here s is a calculated burn-in set of hi/lo-fi simulation pairs
    kd = Dict(kwargs)
    cloud = MakeMFABCCloud(s, epsilons, (1.0, 1.0))
    burnin=length(s)
    if haskey(kd, :N)
        NN = Integer(kd[:N])-burnin # Input specifies the total number of particles (takes precedence over budget), NN here calculates additional particles needed
        for n in 1:NN
            etas = get_eta(cloud, epsilons; method=method, kwargs...)[1]
            push!(cloud, MFABCParticle(mfabc, epsilons, etas)[1])
        end
    elseif haskey(kd, :budget)
        b = Float64(kd[:budget])
        run_cost = cost(s)
        while run_cost<b
            etas = get_eta(cloud, epsilons; method=method, kwargs...)[1]
            pp, c = MFABCParticle(mfabc, epsilons, etas)
            run_cost += c
            push!(cloud, pp)
        end
    end
    return cloud
end

# convert (fixed eta or unknown eta: determined by specifying a burnin or otherwise)

function MakeMFABCCloud(s::BenchmarkCloud, epsilons::Tuple{Float64, Float64}; method::String="mf", kwargs...)
# Continuation probabilities are either inferred from the entire benchmark set (posthoc), or sequentially (after a burn-in period)
    kd = Dict(kwargs)
    if haskey(kd,:burnin)
        n_b = kd[:burnin]::Int64
        return MakeMFABCCloud_(s, epsilons, n_b; method=method, kwargs...)
    else
    # Convert a benchmark cloud into a MFABC cloud using the entire benchmark set's information to get an optimal eta
        etas = get_eta(s, epsilons; method=method, kwargs...)[1]
        return MakeMFABCCloud(s, epsilons, etas)
    end
end
function MakeMFABCCloud_(s::BenchmarkCloud, epsilons::Tuple{Float64, Float64}, burnin_size::Int64; method::String="mf", kwargs...)
# Convert a benchmark cloud into a MFABC cloud sequentially, using only the preceding information to inform eta at each iteration
# This is a very artificial approach, as we are ignoring information from simulations already carried out
    out = MakeMFABCCloud(s[1:burnin_size], epsilons, (1.0,1.0))

    for n in (burnin_size+1):length(s)
        etas = get_eta(out, epsilons; method=method, kwargs...)[1]
        push!(out, MFABCParticle(s[n], epsilons, etas))
    end
    return out
end

########## MFABC Clouds relies on "GET_ETA", calculated as follows:

function sample_properties(s::MFABCCloud, epsilons::Tuple{Float64, Float64}; kwargs...)
    kd = Dict(kwargs)

    if haskey(kd, :F)
        fun = kd[:F]::Function
        F_i = map(pp->fun(pp.p.k), s)
        w_i = [pp.w for pp in s]
        Fbar = sum(F_i .* w_i)/sum(w_i)
        Fweights = (F_i .- Fbar).^2
    else
        Fweights = ones(length(s))
    end

    is_bm = [isa(pp.p, Particle{2}) for pp in s]
    s2 = s[is_bm]
    Fweights2 = Fweights[is_bm]

    rhom = mean([pp.p.dist[1]<epsilons[1] for pp in s])
    rhok = mean([pp.p.dist[1]<epsilons[1] for pp in s2])

    O_lo = [pp.p.dist[1]<epsilons[1] for pp in s2]
    O_hi = [pp.p.dist[2]<epsilons[2] for pp in s2]

    p_tp = (rhom/rhok) * mean(Fweights2.*(O_lo .& O_hi))
    p_fp = (rhom/rhok) * mean(Fweights2.*(O_lo .& .~O_hi))
    p_fn = ((1-rhom)/(1-rhok)) * mean(Fweights2.*(.~O_lo .& O_hi))
    ct = mean([pp.p.cost[1] for pp in s])
    c_p = (rhom/rhok) * mean([pp.p.cost[2] for pp in s2] .* O_lo)
    c_n = ((1-rhom)/(1-rhok)) * mean([pp.p.cost[2] for pp in s2] .* .~O_lo)

    return p_tp, p_fp, p_fn, ct, c_p, c_n
end

function sample_properties(s::BenchmarkCloud, epsilons::Tuple{Float64, Float64}; kwargs...)
    kd = Dict(kwargs)
    if haskey(kd, :F)
        fun = kd[:F]::Function
        F_i = map(p->fun(p.k), s)
        w_i = [p.dist[2]<epsilons[2] for p in s]
        Fbar = sum(F_i .* w_i)/sum(w_i)
        Fweights = (F_i .- Fbar).^2
    else
        Fweights = ones(length(s))
    end

    O_lo = [p.dist[1]<epsilons[1] for p in s]
    O_hi = [p.dist[2]<epsilons[2] for p in s]

    p_tp = mean(Fweights.*(O_lo .& O_hi))
    p_fp = mean(Fweights.*(O_lo .& .~O_hi))
    p_fn = mean(Fweights.*(.~O_lo .& O_hi))
    ct = mean([p.cost[1] for p in s])
    c_p = mean([p.cost[2] for p in s] .* O_lo)
    c_n = mean([p.cost[2] for p in s] .* .~O_lo)

    return p_tp, p_fp, p_fn, ct, c_p, c_n
end

function phi(eta, p_tp::Float64, p_fp::Float64, p_fn::Float64, ct::Float64, c_p::Float64, c_n::Float64)
    return (p_tp - p_fp + (p_fp/eta[1]) + (p_fn/eta[2])) * (ct + c_p*eta[1] + c_n*eta[2])
end
function phi(eta, s::Cloud, epsilons::Tuple{Float64, Float64}; kwargs...)
    return phi(eta, sample_properties(s, epsilons; kwargs...)...)
end

function get_eta(p_tp::Float64, p_fp::Float64, p_fn::Float64, ct::Float64, c_p::Float64, c_n::Float64; method::String="mf")
    Rp = p_fp/(c_p/ct)
    Rn = p_fn/(c_n/ct)
    R0 = p_tp - p_fp
    if R0<=0
        return (1.0,1.0)
    end

    eta_1 = sqrt(Rp/R0)
    eta_2 = sqrt(Rn/R0)

    eta_1_bar = minimum([1.0, eta_1 / sqrt((1+p_fn/R0)/(1+c_n/ct))])
    eta_2_bar = minimum([1.0, eta_2 / sqrt((1+p_fp/R0)/(1+c_p/ct))])

    if method=="mf"
        if (eta_1<=1.0)&(eta_2<=1.0)
            return (eta_1,eta_2)
        else
            phi_1 = phi((eta_1_bar, 1.0), p_tp, p_fp, p_fn, ct, c_p, c_n)
            phi_2 = phi((1.0, eta_2_bar), p_tp, p_fp, p_fn, ct, c_p, c_n)
            if phi_2 < phi_1
                return (1.0, eta_2_bar)
            else
                return (eta_1_bar, 1.0)
            end
        end
    elseif method=="er"
        return (1.0, eta_2_bar)
    elseif method=="ed"
        eta = minimum([1.0, sqrt((ct*(p_fp+p_fn))/(R0*(c_p+c_n)))])
        return (eta, eta)
    else
        return (1.0, 1.0)
    end
end

function get_eta(s::Cloud, epsilons::Tuple{Float64, Float64}; method::String="mf", lower_eta::Float64=0.01, kwargs...)
    # Keyword "method" can be chosen from "mf" (multifidelity), "er" (early rejection), "ed" (early decision)
    # Any other keyword will give the ABC approach with no early rejection or acceptance at all.

    sp = sample_properties(s, epsilons; kwargs...)
    eta = get_eta(sp..., method=method)
    eta = max.(lower_eta, eta)
    return eta, phi(eta, sp...)
end

######## Read/write (benchmark) clouds

using DelimitedFiles

function write_cloud(c::Array, name::String)
    open(name*".bm", "w") do io
        writedlm(io, c)
    end
    return nothing
end
function write_cloud(c::Array{<:Particle,1}, name::String)
    mkpath(name)
    for fieldname in fieldnames(Particle)
        cloud_field = [getfield(p, fieldname) for p in c]
        write_cloud(cloud_field, name*"/"*string(fieldname))
    end
    return nothing
end
function write_cloud(c::MFABCCloud, name::String)
    for fieldname in fieldnames(MFABCParticle)
        cloud_field = [getfield(mfp, fieldname) for mfp in c]
        write_cloud(cloud_field, name*"/"*string(fieldname))
    end
    N = length.(c)
    write_cloud(N, name*"/p/N")
    return nothing
end

function MakeBenchmarkCloud(indir::String)::BenchmarkCloud
# Read a previously simulated Benchmark Cloud
    k = readdlm(indir*"/k.bm")
    dist = readdlm(indir*"/dist.bm")
    cost = readdlm(indir*"/cost.bm")
    n = size(k,1)

    bm = BenchmarkCloud(undef, n)
    for i in 1:n
        bm[i] = Particle{2}(Parameters(k[i,:]), NTuple{2,Float64}(dist[i,:]), NTuple{2,Float64}(cost[i,:]))
    end
    return bm
end

function MakeMFABCCloud(indir::String)::MFABCCloud
    N = Integer.(readdlm(indir*"/p/N.bm"))
    k = readdlm(indir*"/p/k.bm")
    dist = readdlm(indir*"/p/dist.bm")
    cost = readdlm(indir*"/p/cost.bm")
    eta = readdlm(indir*"/eta.bm")
    w = readdlm(indir*"/w.bm")

    n = size(k,1)
    mf = MFABCCloud(undef, n)
    for i in 1:n
        pp = Particle{N[i]}(Parameters(k[i,:]), NTuple{N[i],Float64}(dist[i,1:N[i]]), NTuple{N[i],Float64}(cost[i,1:N[i]]))
        mf[i] = MFABCParticle(pp, eta[i], w[i])
    end
    return mf
end

max_stars_repo_name: tpprescott/mf-abc

max_issues_repo_name: tpprescott/mf-abc

max_forks_repo_name: tpprescott/mf-abc


text: draw(widget, area, buffer) = nothing


text:
using Test
using VisualDL
using PyCall
@pyimport visualdl as vdl

train_logger = VisualDLLogger("tmp", 1, "train")
test_logger = as_mode(train_logger, "test")

for i in 1:3
    with_logger(train_logger) do
        @log_scalar s0=(i,i)
    end

    with_logger(test_logger) do
        @log_scalar s1=(i, -i)
    end
end

set_caption(train_logger, :s0, "This is caption of s0")
set_caption(test_logger, :s1, "This is caption of s1")

save(train_logger)
save(test_logger)

@testset "Test VisualDLLogger" begin
    reader = vdl.LogReader("tmp")

    @test reader[:modes]() == ["train", "test"]

    reader[:mode]("train")
    @test reader[:tags]("scalar") == ["s0"]
    s = reader[:scalar]("s0")
    @test s[:caption]() == "This is caption of s0"
    @test s[:ids]() == [1, 2, 3]
    @test s[:records]() == [1, 2, 3]

    reader[:mode]("test")
    @test reader[:tags]("scalar") == ["s1"]
    s = reader[:scalar]("s1")
    @test s[:ids]() == [1, 2, 3]
    @test s[:records]() == [-1, -2, -3]
end

rm("tmp"; recursive=true)


text:
struct QBosonicBasis <: AbstractQBasis
  	Dim::Int
	N::Int
	Np::Int
	OccRep::Vector{Vector{Int}}
end

function QBosonicBasis(N::Int, Np::Int)
	I = 1
	n = zeros(Int,N)
	n[1] = Np

	vectors = Vector{Vector{Int}}([])
	while true
		#@show n
		push!(vectors,n[:])
		I += 1
		if n[1] > 0
			n[1] -= 1
			n[2] += 1
		elseif n[1] == 0
			j = findfirst(!iszero, n)

			if j == N break end

			n[1] = n[j] - 1
			n[j] = 0
			n[j+1] += 1

		end
		#@show n
	end

	return QBosonicBasis(I-1, N, Np, vectors)
end

function bosonicHashfunc(Np::Int, col::Vector{Int})
	hashkey = 0
	particlei = 1
	p_i = findfirst(!iszero,col)
	col = col[p_i:end]
	while true

		#@show col

		if col[1] > 0
			col[1] -= 1
		else
 			orbital = findfirst(!iszero,col)
			p_i += orbital-1
			col = col[orbital:end]
			col[1] -= 1
		end

		hashkey += binomial((p_i-1) + particlei -1 , particlei)

		if particlei == Np break end
 		particlei += 1
	end
	return  hashkey
end

function addParticleBosonic(singleCoef::AbstractVector{<:Number}, preState::QState)
 	N = preState.QBasis.N
	Np = preState.QBasis.Np

	N >= Np+1 || throw(DomainError("you cannot creat more particles $(Np+1) than sites $N"))

	newQBasis = QBosonicBasis(N,Np+1)
	newCoef = zeros(ComplexF64,newQBasis.Dim)
	preQBasis = preState.QBasis

	for i in 1:N
	 	for (colIdx, colOccRep) in enumerate(preQBasis.OccRep)
			#@show colOccRep
	 		colOccRep[i] += 1

			# actualy this bosonicHashfunc gives directly idx in
			# OccRep vector so no dict needed
			rowIdx = bosonicHashfunc(Np+1, colOccRep)+1
			#@show i,rowIdx
			newCoef[rowIdx] +=  singleCoef[i]*preState.Coef[colIdx]*sqrt(colOccRep[i])

			#@show newCoef[rowIdx]
			colOccRep[i] -= 1
		end
		#@show newCoef
	end

	return QState(newCoef, newQBasis)
end

function constructBosonicState(blochWaveFunc::AbstractVector{<:Any})
	#@show size(blochWaveFunc[1,:,1])
	State = QState(blochWaveFunc[1][:,1], QBosonicBasis(size(blochWaveFunc[1])[1],1))
	for i in 2:size(blochWaveFunc)[1]
		State = addParticleBosonic(blochWaveFunc[i][:,1], State)
	end

	return State
end

max_stars_repo_name: Gawatz/exactQmodule

max_issues_repo_name: Gawatz/exactQmodule

max_forks_repo_name: Gawatz/exactQmodule


text:
# Autogenerated wrapper script for Libiconv_jll for i686-linux-musl
export libcharset, libiconv

JLLWrappers.@generate_wrapper_header("Libiconv")
JLLWrappers.@declare_library_product(libcharset, "libcharset.so.1")
JLLWrappers.@declare_library_product(libiconv, "libiconv.so.2")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcharset,
        "lib/libcharset.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libiconv,
        "lib/libiconv.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Plots
using LinearAlgebra
using Distributions

xC = [2, 1] # Center of data (mean)
sig = [2, 0.5] # Principal axes

theta = pi/3 # Rotate cloud by pi/3
R = [cos(theta) -sin(theta); sin(theta) cos(theta)] # Rotation matrix

nPoints = 10000 # Create 10,000 points
X = R*diagm(sig)*rand(Normal(), 2, nPoints) + diagm(xC)*ones(2, nPoints);

p1 = scatter(X[1,:], X[2,:]) # Plot cloud of noisy data

# f_ch01_ex03_1b

Xavg = [sum(X[1,:]); sum(X[2,:])]/nPoints # Compute mean
B = X - Xavg*ones(1,nPoints) # Mean-subtracted Data
U, S, V = svd(B/sqrt(nPoints)) # Find principal components (SVD)

theta = 2pi*(0:0.01:1)
Xstd = hcat([U*[S[1]*cos(t); S[2]*sin(t)] for t in theta]...) # 1-std confidence interval
p2 = plot(Xstd[1,:] .+ Xavg[1], Xstd[2,:] .+ Xavg[2])
plot!(2*Xstd[1,:] .+ Xavg[1], 2*Xstd[2,:] .+ Xavg[2])
plot!(3*Xstd[1,:] .+ Xavg[1], 3*Xstd[2,:] .+ Xavg[2])

# Plot principal components
# p3 = scatter([Tuple(Xavg + U[:,1]*S[1]), Tuple(Xavg + U[:,2]*S[2])])


text:
#!/usr/bin/env julia
"""
Testing code for other scripts

"""
using VariantCall

chr = ARGS[1]
inputf = ARGS[2]

line = VariantCall.read(inputf)
record = VariantCall.records(inputf)
println(join(record,"\n"))
print(join(record[1:15]), "\n")

#inf = VariantCall.info(record)
#aa = VariantCall.infoAA(record)
thouG = VariantCall.E1000G(record)
#println(info)

#from, to = chrange(record, chr)
#extract = VariantCall.fetch(record, chr, 17104729, 18078510)

#alterations = VariantCall.uniqalts(record)
#for ln in alterations#[1:30]
#	println(ln)
#end
the stupid way of programming


text:
@time using FastGroupBy
@time sumby(id6,v1)
@time addprocs()

@everywhere using FastGroupBy
@everywhere using SplitApplyCombine
@time psumby(id6,v1)

@time as = SharedArray(id6);
@time bs = SharedArray(v1);

@time psumby(as,bs)

@time pgroupreduce(x -> x[1], x->x[2], (x,y)-> x+y, dict_add_reduce, as, bs); #29
@time pgroupreduce(x -> x[1], x->x[2], (x,y)-> x+y, dict_add_reduce, as, bs)
@time psumby(id6, v1); #26

@time res = pgroupreduce(x -> x[1], x->(x[2],1), (x,y)-> (x[1]+y[1], x[2]+y[2]), dict_mean_reduce, as, bs); #29

@time countmap(id6,weights(v1))
@time sumby(id6,v1)
@time groupreduce(x->x[1], x->x[2], (x,y)->x+y, zip(id6,v1))


text:
using Test
using ClimaCore.Geometry, ClimaCore.DataLayouts
using LinearAlgebra, StaticArrays

@testset "AxisTensors" begin
    x = Geometry.Covariant12Vector(1.0, 2.0)
    y = Geometry.Contravariant12Vector(1.0, 4.0)

    M = Geometry.Axis2Tensor(
        (Geometry.Cartesian12Axis(), Geometry.Covariant12Axis()),
        [1.0 0.0; 0.5 2.0],
    )

    @test dot(x, y) == x' * y == 9.0
    @test dot(y, x) == y' * x == 9.0

    @test x == x
    @test x != Geometry.components(x)
    @test x != Geometry.Contravariant12Vector(Geometry.components(x)...)

    @test x[1] == 1.0
    @test y[2] == 4.0
    @test M[2] == 0.5
    @test M[2, 1] == 0.5
    @test M[:, 1] == Geometry.Cartesian12Vector(1.0, 0.5)
    @test M[1, :] == Geometry.Covariant12Vector(1.0, 0.0)

    @test x * y' ==
          x ⊗ y ==
          Geometry.AxisTensor(
              (axes(x, 1), axes(y, 1)),
              Geometry.components(x) * Geometry.components(y)',
          )

    @test Geometry.components(M * inv(M)) == @SMatrix [1.0 0.0; 0.0 1.0]
    @test Geometry.components(inv(M) * M) == @SMatrix [1.0 0.0; 0.0 1.0]

    @test x ⊗ 3 == Geometry.Covariant12Vector(3.0, 6.0)
    @test x ⊗ (1, (a = 2, b = 3)) == (
        Geometry.Covariant12Vector(1.0, 2.0),
        (
            a = Geometry.Covariant12Vector(2.0, 4.0),
            b = Geometry.Covariant12Vector(3.0, 6.0),
        ),
    )

    @test Geometry.components(M * inv(M)) == @SMatrix [1.0 0.0; 0.0 1.0]
    @test Geometry.components(inv(M) * M) == @SMatrix [1.0 0.0; 0.0 1.0]

    @test M * y == Geometry.Cartesian12Vector(1.0, 8.5)
    @test M \ Geometry.Cartesian12Vector(1.0, 8.5) == y

    @test_throws DimensionMismatch dot(x, x)
    @test_throws DimensionMismatch M * x
    @test_throws DimensionMismatch M \ x

    @test DataLayouts.typesize(Float64, typeof(x)) == 2
end

@testset "transform" begin
    @test Geometry.transform(
        Geometry.Covariant12Axis(),
        Geometry.Covariant1Vector(2.0),
    ) == Geometry.Covariant12Vector(2.0, 0.0)
    @test Geometry.transform(
        Geometry.Covariant12Axis(),
        Geometry.Covariant13Vector(2.0, 0.0),
    ) == Geometry.Covariant12Vector(2.0, 0.0)
    @test_throws InexactError Geometry.transform(
        Geometry.Covariant12Axis(),
        Geometry.Covariant13Vector(2.0, 2.0),
    )

    @test Geometry.transform(
        Geometry.Covariant12Axis(),
        Geometry.Covariant1Vector(2.0) * Geometry.Cartesian1Vector(1.0)',
    ) == Geometry.Covariant12Vector(2.0, 0.0) * Geometry.Cartesian1Vector(1.0)'
    @test Geometry.transform(
        Geometry.Covariant12Axis(),
        Geometry.Covariant13Vector(2.0, 0.0) * Geometry.Cartesian1Vector(1.0)',
    ) == Geometry.Covariant12Vector(2.0, 0.0) * Geometry.Cartesian1Vector(1.0)'
    @test_throws InexactError Geometry.transform(
        Geometry.Covariant12Axis(),
        Geometry.Covariant13Vector(2.0, 2.0) * Geometry.Cartesian1Vector(1.0)',
    )
end

@testset "project" begin
    @test Geometry.project(
        Geometry.Covariant12Axis(),
        Geometry.Covariant1Vector(2.0),
    ) == Geometry.Covariant12Vector(2.0, 0.0)

    @test Geometry.project(
        Geometry.Covariant12Axis(),
        Geometry.Covariant12Vector(2.0, 2.0),
    ) == Geometry.Covariant12Vector(2.0, 2.0)

    @test Geometry.project(
        Geometry.Covariant12Axis(),
        Geometry.Covariant123Vector(2.0, 2.0, 0.0),
    ) == Geometry.Covariant12Vector(2.0, 2.0)

    @test Geometry.project(
        Geometry.Covariant12Axis(),
        Geometry.Covariant123Vector(2.0, 2.0, 1.0),
    ) == Geometry.Covariant12Vector(2.0, 2.0)

    @test Geometry.project(
        Geometry.Covariant12Axis(),
        Geometry.Covariant13Vector(2.0, 0.0),
    ) == Geometry.Covariant12Vector(2.0, 0.0)
    @test Geometry.project(
        Geometry.Covariant12Axis(),
        Geometry.Covariant13Vector(2.0, 2.0),
    ) == Geometry.Covariant12Vector(2.0, 0.0)

    @test Geometry.project(
        Geometry.Covariant12Axis(),
        Geometry.Covariant1Vector(2.0) * Geometry.Cartesian1Vector(1.0)',
    ) == Geometry.Covariant12Vector(2.0, 0.0) * Geometry.Cartesian1Vector(1.0)'
    @test Geometry.project(
        Geometry.Covariant12Axis(),
        Geometry.Covariant13Vector(2.0, 0.0) * Geometry.Cartesian1Vector(1.0)',
    ) == Geometry.Covariant12Vector(2.0, 0.0) * Geometry.Cartesian1Vector(1.0)'
    @test Geometry.project(
        Geometry.Covariant12Axis(),
        Geometry.Covariant13Vector(2.0, 2.0) * Geometry.Cartesian1Vector(1.0)',
    ) == Geometry.Covariant12Vector(2.0, 0.0) * Geometry.Cartesian1Vector(1.0)'
end

@testset "cross product" begin
    M = @SMatrix [
        4.0 1.0
        0.5 2.0
    ]
    J = det(M)
    local_geom = Geometry.LocalGeometry(
        Geometry.XYPoint(0.0, 0.0),
        J,
        J,
        Geometry.Axis2Tensor(
            (Geometry.UVAxis(), Geometry.Covariant12Axis()),
            M,
        ),
    )

    u = Geometry.UVVector(1.0, 2.0)
    v = Geometry.WVector(3.0)
    @test u × v == -v × u == Geometry.UVVector(6.0, -3.0)
    uⁱ = Geometry.ContravariantVector(u, local_geom)
    vⁱ = Geometry.ContravariantVector(v, local_geom)
    @test Geometry.UVVector(Geometry._cross(uⁱ, vⁱ, local_geom), local_geom) ==
          Geometry.UVVector(6.0, -3.0)
end