text:
module TesSpaces

using Test
using TreeParzen

space = Dict(
    :a => HP.Choice(:a, [:a, :b, :c]),
    :b => HP.Choice(:b, ["a", "b", "c"])
)

# just asking, isn't good enough
trial = ask(space)
sample = trial.hyperparams

@test sample[:a] in (:a, :b, :c)
@test sample[:b] in ("a", "b", "c")

# gotta ask both ways
trials = TreeParzen.Trials.Trial[]
tell!(trials, trial, 2.)
trial = ask(space)
tell!(trials, trial, 2.6)

tp_trial = ask(space, trials, TreeParzen.Config())
tp_sample = tp_trial.hyperparams

# same holds true but we have to make sure we can sample from strings literals via TP asks
@test tp_sample[:a] in (:a, :b, :c)
@test tp_sample[:b] in ("a", "b", "c")

end
true

ext: jl


text:
using LeastSquaresEstimation
using Plots

function main()
    # x = [k₁, k₂, k₃, k₄, k₅, λ₁, λ₂, λ₃, ω₁, ω₂, ω₃, δ₁, δ₂, δ₃]ᵀ

    # Define basis function vector
    θ(x,t) = x[1]*exp(x[6]*t)*cos(x[9]*t + x[12]) +
             x[2]*exp(x[7]*t)*cos(x[10]*t + x[13]) +
             x[3]*exp(x[8]*t)*cos(x[11]*t + x[14]) + x[4]
    ϕ(x,t) = x[1]*exp(x[6]*t)*sin(x[9]*t + x[12]) +
             x[2]*exp(x[7]*t)*sin(x[10]*t + x[13]) +
             x[3]*exp(x[8]*t)*sin(x[11]*t + x[14]) + x[5]
    f = [θ, ϕ]

    # Instantiate Linear Least Squares Estimator
    xhat0 = [0.5, 0.25, 0.125, 0.0, 0.0, -0.15, 0.06,
            -0.03, 0.26, 0.55, 0.95, 0.01, 0.01, 0.01]
    nllse = NonLinearLeastSquares(f, xhat0)

    # Generate Measurements
    ft(x,t) = [f[1](x,t), f[2](x,t)] .+ 0.0005.*randn(2)
    m    = 200
    p    = length(f)
    xt   = [0.2, 0.1, 0.05, 0.0001, 0.0001, -0.1, -0.05,
           -0.025, 0.25, 0.5, 1.0, 0.0, 0.0, 0.0]

    ml   = Int(m/p)
    ts = zeros(ml)
    ys = zeros(m)
    @inbounds for i in 1:ml
        if i != 1
            ts[i] = ts[i-1] + 0.1
        end
        idxs = ((i-1)*2 + 1):((i-1)*2 + p)
        ys[idxs] = ft(xt, ts[i])
    end

    # Feed estimator batch of measurements
    W = zeros(m,m)
    @inbounds for i in 1:m
        W[i,i] = 1.0
    end
    FeedMeasurementBatch!(nllse, ys, W, ts)

    # Compute estimate
    ComputeEstimate!(nllse)

    # Compute estimated measurements and plot
    yhat = zeros(m)
    @inbounds for i in 1:ml
        yhat[(i-1)*p + 1:(i-1)*p + p] .= GetMeasurementEstimate(nllse, ts[i])
    end

    # Display
    println("Estimated parameters:")
    display(GetEstimate(nllse))

    # Plot generated and estimated measurements
    p1 = plot(ts, ys[1:2:end]; label = "Simulated θ")
    plot!(p1, ts, yhat[1:2:end]; label = "Estimated θ")

    p2 = plot(ts, ys[2:2:end]; label = "Simulated ϕ")
    plot!(p2, ts, yhat[2:2:end]; label = "Estimated ϕ")

    p3 = plot(ts, ys[1:2:end] .- yhat[1:2:end]; label = "θ Estimation Error")
    p4 = plot(ts, ys[2:2:end] .- yhat[2:2:end]; label = "ϕ Estimation Error")

    plot(p1,p2,p3,p4,layout = (4,1))
end

main()


text:
"""
clip_bbox(x::AbstractGeoArray, b::bbox)
clip_bbox(x::AbstractGeoArray, y::AbstractGeoArray)
"""
function st_clip(ga::AbstractGeoArray, b::bbox)
    x, y = st_dim(ga)
    I_x = findall((x .<= b.xmax) .& (x .>= b.xmin))
    I_y = findall((y .<= b.ymax) .& (y .>= b.ymin))
    ga[I_x, I_y]
end

function st_clip(x::AbstractGeoArray, y::AbstractGeoArray)
    st_clip(x, st_bbox(y))
end

export st_clip


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

@test convert(Tuple, (1,2)) == (1,2)
@testset "indexing" begin
    @test length(()) == 0
    @test length((1,)) === 1
    @test length((1,2)) === 2

    @test_throws ArgumentError Base.front(())
    @test_throws ArgumentError first(())

    @test endof(()) === 0
    @test endof((1,)) === 1
    @test endof((1,2)) === 2

    @test size((), 1) === 0
    @test size((1,), 1) === 1
    @test size((1,2), 1) === 2
end
# @test_throws ArgumentError size((), 2)
# @test_throws ArgumentError size((1,), 2)
# @test_throws ArgumentError size((1,2), 2)

@test getindex((1,), 1) === 1
@test getindex((1,2), 2) === 2
@test_throws BoundsError getindex((), 1)
@test_throws BoundsError getindex((1,2), 0)
@test_throws BoundsError getindex((1,2), -1)

@test getindex((1,), 1.0) === 1
@test getindex((1,2), 2.0) === 2
@test_throws BoundsError getindex((), 1.0)
@test_throws BoundsError getindex((1,2), 0.0)
@test_throws BoundsError getindex((1,2), -1.0)

@test getindex((5,6,7,8), [1,2,3]) === (5,6,7)
@test_throws BoundsError getindex((1,2), [3,4])

@test getindex((5,6,7,8), [true, false, false, true]) === (5,8)
@test_throws BoundsError getindex((5,6,7,8), [true, false, false, true, true])

@test getindex((5,6,7,8), []) === ()

## filling to specified length
@test @inferred(Base.fill_to_length((1,2,3), -1, Val{5})) == (1,2,3,-1,-1)
@test_throws ArgumentError Base.fill_to_length((1,2,3), -1, Val{2})

## iterating ##
@test start((1,2,3)) === 1

@test done((), 1)
@test !done((1,2,3), 3)
@test done((1,2,3), 4)

@test next((5,6,7), 1) === (5, 2)
@test next((5,6,7), 3) === (7, 4)
@test_throws BoundsError next((5,6,7), 0)
@test_throws BoundsError next((), 1)

@test collect(eachindex((2,5,"foo"))) == collect(1:3)
@test collect(eachindex((2,5,"foo"), (1,2,5,7))) == collect(1:4)

## eltype ##
@test eltype((1,2,3)) === Int
@test eltype((1.0,2.0,3.0)) <: AbstractFloat
@test eltype((true, false)) === Bool
@test eltype((1,2.0, false)) === Any
@test eltype(()) === Union{}

begin
    local foo
    ## mapping ##
    foo() = 2
    foo(x) = 2x
    foo(x, y) = x + y
    foo(x, y, z) = x + y + z
    longtuple = ntuple(identity, 20)

    # 1 argument
    @test map(foo, ()) === ()
    @test map(foo, (1,)) === (2,)
    @test map(foo, (1,2)) === (2,4)
    @test map(foo, (1,2,3,4)) === (2,4,6,8)
    @test map(foo, longtuple) === ntuple(i->2i,20)

    # 2 arguments
    @test map(foo, (), ()) === ()
    @test map(foo, (1,), (1,)) === (2,)
    @test map(foo, (1,2), (1,2)) === (2,4)
    @test map(foo, (1,2,3,4), (1,2,3,4)) === (2,4,6,8)
    @test map(foo, longtuple, longtuple) === ntuple(i->2i,20)

    # n arguments
    @test map(foo, (), (), ()) === ()
    @test map(foo, (), (1,2,3), (1,2,3)) === ()
    @test map(foo, (1,), (1,), (1,)) === (3,)
    @test map(foo, (1,2), (1,2), (1,2)) === (3,6)
    @test map(foo, (1,2,3,4), (1,2,3,4), (1,2,3,4)) === (3,6,9,12)
    @test map(foo, longtuple, longtuple, longtuple) === ntuple(i->3i,20)
end

## comparison ##
@test isequal((), ())
@test isequal((1,2,3), (1,2,3))
@test !isequal((1,2,3), (1,2,4))
@test !isequal((1,2,3), (1,2))

@test ==((), ())
@test ==((1,2,3), (1,2,3))
@test !==((1,2,3), (1,2,4))
@test !==((1,2,3), (1,2))

@test isless((1,2), (1,3))
@test isless((1,), (1,2))
@test !isless((1,2), (1,2))
@test !isless((2,1), (1,2))

## functions ##
@test isempty(())
@test !isempty((1,))

@test reverse(()) === ()
@test reverse((1,2,3)) === (3,2,1)

## specialized reduction ##
@test sum((1,2,3)) === 6

@test prod(()) === 1
@test prod((1,2,3)) === 6

@test all(()) === true
@test all((false,)) === false
@test all((true,)) === true
@test all((true, true)) === true
@test all((true, false)) === false
@test all((false, false)) === false

@test any(()) === false
@test any((true,)) === true
@test any((false,)) === false
@test any((true, true)) === true
@test any((true, false)) === true
@test any((false, false)) === false
@test any((false,false,false)) === false
@test any((true,false,false)) === true
@test any((false,true,false)) === true
@test any((false,false,true)) === true
@test any((true,true,false)) === true
@test any((true,false,true)) === true
@test any((true,true,false)) === true
@test any((true,true,true)) === true

@test @inferred(ntuple(abs2, Val{0})) == ()
@test @inferred(ntuple(abs2, Val{2})) == (1, 4)
@test @inferred(ntuple(abs2, Val{3})) == (1, 4, 9)
@test @inferred(ntuple(abs2, Val{4})) == (1, 4, 9, 16)
@test @inferred(ntuple(abs2, Val{5})) == (1, 4, 9, 16, 25)
@test @inferred(ntuple(abs2, Val{6})) == (1, 4, 9, 16, 25, 36)

# issue #12854
@test_throws TypeError ntuple(identity, Val{1:2})

for n = 0:20
    t = ntuple(identity, n)
    @test length(t) == n
    for i = 1:n
        @test t[i] == i
    end
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
#   Unit examples of the FourierAnalysis Package for julia language
#   v 0.0.1 - last update 24th of September 2019
#
#   MIT License
#   Copyright (c) 2019, Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This example show how to compute coherence
#   and how to extract inforation from them.

#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #
#                                                                             #
#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #

using FourierAnalysis, FFTW, LinearAlgebra, Statistics, Plots

# add module for reading the two EEG text files to be used ater
push!(LOAD_PATH, @__DIR__)
using IOtxt

# Get EEG file names with complete path (they have extension .txt)
S=getFilesInDir(@__DIR__; ext=(".txt",))

# read the two EEG data files and put them in a Matrix object
X1=readEEG(S[1])
X2=readEEG(S[2])

#############################################################
# test coherence using sinusoids with same frequency and different phase.
using FourierAnalysis, FFTW, LinearAlgebra, Statistics, Plots
t=32
sr=32
f=2
noise=1000
tapering=rectangular
v=sinusoidal(10, f, sr, t*4, 0).+[randn()/noise for i=1:t*4]
w=sinusoidal(10, f, sr, t*4, π/4).+[randn()/noise for i=1:t*4]
y=sinusoidal(10, f, sr, t*4, π/2).+[randn()/noise for i=1:t*4]
z=sinusoidal(10, f, sr, t*4, π).+[randn()/noise for i=1:t*4]
X=[v w y z]
# as compared to v, w is out-of-phase, y is in between and z is at phase opposition
plot(v)
plot!(w)
plot!(y)
plot!(z)

𝘾=coherence(X, sr, t; tapering=tapering, tril=true)
𝘾=coherence(X, sr, t; tapering=tapering, tril=false)
# get all five kinds of coherences
𝘾₁, 𝘾₂, 𝘾₃, 𝘾₄, 𝘾₅=coherence(X, sr, t; tapering=tapering, tril=true, allkinds=true)
𝘾₁, 𝘾₂, 𝘾₃, 𝘾₄, 𝘾₅=coherence(X, sr, t; tapering=tapering, tril=false, allkinds=true)
# check result against cross-spectra
𝐒=crossSpectra(X, sr, t; tapering=tapering, tril=true)
# compute coherence from a CrossSpectra object (lower triangles only)
𝘾=coherence(𝐒)

𝐒=crossSpectra(X, sr, t; tapering=tapering, tril=false)
# compute coherence from a CrossSpectra object (full matrices)
𝘾=coherence(𝐒)

# compute all five kinds of coherence from a CrossSpectra object
𝘾₁, 𝘾₂, 𝘾₃, 𝘾₄, 𝘾₅=coherence(𝐒, allkinds=true)

# average full coherence in range 4:8Hz
Y=mean(𝘾, (4, 8))

# full coherence matrices in range 4:8Hz
E=extract(𝘾, (4, 8))

# average coherence in 1Hz-bands
Z=bands(𝘾, 1)

# smooth coherence a-posteriori
𝘾2=smooth(blackmanSmoother, 𝘾)
C=𝘾.y
D=𝘾2.y

# or get directly smoothed coherence
𝘾3=coherence(X, sr, t; tapering=tapering, smoothing=hannSmoother, tril=true)

############################################
# coherence of several data matrices at once

# Compute the coherence altogether
𝓒=coherence(𝐗, sr, t; tapering=slepians(sr, t))

# get all five kinds of coherences
𝓒₁, 𝓒₂, 𝓒₃, 𝓒₄, 𝓒₅=coherence(𝐗, sr, t; tapering=slepians(sr, t), allkinds=true)


text:
using DistStat, Random, LinearAlgebra

mutable struct COXUpdate
    maxiter::Int
    step::Int
    verbose::Bool
    tol::Real
    function COXUpdate(; maxiter::Int=100, step::Int=10, verbose::Bool=false, tol::Real=1e-10)
        maxiter > 0 || throw(ArgumentError("maxiter must be greater than 0."))
        tol > 0 || throw(ArgumentError("tol must be positive."))
        new(maxiter, step, verbose, tol)
    end
end

mutable struct COXVariables{T, A}
    # COX variables corresponding to an m × [n] data matrix X (external). m: number of subjects, n: number of predictors
    m::Int # rows, number of subjects
    n::Int # cols, number of predictors
    β::MPIVector{T,A} # [n]
    β_prev::MPIVector{T,A} # r x [n]
    δ::A # indicator for right censoring (0 if censored)
    λ::T # regularization parameter
    t::AbstractVector{T} # a vector containing timestamps, should be in descending order.
    σ::T # step size. 1/(2 * opnorm(X)^2) for guaranteed convergence
    π_ind::MPIMatrix{T,A}
    tmp_n::MPIVector{T,A}
    tmp_1m::MPIMatrix{T,A}
    tmp_mm::MPIMatrix{T,A}
    tmp_m_local1::A
    tmp_m_local2::A
    eval_obj::Bool
    function COXVariables(X::MPIMatrix{T,A}, δ::A, λ::AbstractFloat,
                            t::A;
                            σ::T=convert(T, 1/(2*opnorm(X; verbose=true)^2)),
                            eval_obj=false) where {T,A}
        m, n = size(X)
        β = MPIVector{T,A}(undef, n)
        β_prev = MPIVector{T,A}(undef, n)
        fill!(β, zero(T))
        fill!(β_prev, zero(T))

        δ = convert(Vector{T}, δ)

        π_ind = MPIMatrix{T,A}(undef, m, m)
        t_dist = distribute(reshape(t, 1, :))
        fill!(π_ind, one(T))
        π_ind .= ((π_ind .* t_dist) .- t) .<= 0

        tmp_n = MPIVector{T,A}(undef, n)
        tmp_1m = MPIMatrix{T,A}(undef, 1, m)
        tmp_mm = MPIMatrix{T,A}(undef, m, m)
        tmp_m_local1 = A{T}(undef, m)
        tmp_m_local2 = A{T}(undef, m)

        new{T,A}(m, n, β, β_prev, δ, λ, t, σ, π_ind, tmp_n, tmp_1m, tmp_mm, tmp_m_local1, tmp_m_local2, eval_obj)
    end
end

function reset!(v::COXVariables{T,A}; seed=nothing) where {T,A}
    fill!(v.β, zero(T))
    fill!(v.β_prev, zero(T))
end

function soft_threshold(x::T, λ::T) ::T where T <: AbstractFloat
    x > λ && return (x - λ)
    x < -λ && return (x + λ)
    return zero(T)
end

function update!(X::MPIArray, u::COXUpdate, v::COXVariables{T,A}) where {T,A}
    mul!(v.tmp_m_local1, X, v.β) # {m} = {m x [n]} * {[n]}.
    v.tmp_m_local1 .= exp.(v.tmp_m_local1) # w
    cumsum!(v.tmp_m_local2, v.tmp_m_local1) # W. TODO: deal with ties.
    v.tmp_1m .= distribute(reshape(v.tmp_m_local2, 1, :)) # W_dist: distribute W.
    v.tmp_mm .= v.π_ind .* v.tmp_m_local1 ./ v.tmp_1m # (π_ind .* w) ./ W_dist. computation order is determined for CuArray safety.
    pd = mul!(v.tmp_m_local1, v.tmp_mm, v.δ) # {m} = {m x [m]} * {m}.
    v.tmp_m_local2 .= v.δ .- pd # {m}.
    grad = mul!(v.tmp_n, transpose(X), v.tmp_m_local2) # {[n]} = {[n] x m} * {m}.
    v.β .= soft_threshold.(v.β .+ v.σ .* grad , v.λ) # {[n]}.
end

function get_objective!(X::MPIArray, u::COXUpdate, v::COXVariables{T,A}) where {T,A}
    v.tmp_n .= (v.β .== 0)
    sparsity = sum(v.tmp_n)/size(v.tmp_n, 1)
    if v.eval_obj
        v.tmp_m_local1 .= exp.(mul!(v.tmp_m_local1, X, v.β))
        cumsum!(v.tmp_m_local1, v.tmp_m_local1)
        obj = dot(v.δ, mul!(v.tmp_m_local2, X, v.β) .- log.(v.tmp_m_local1)) .- v.λ .* sum(abs.(v.β)) # TODO: deal with ties.

        return false, (obj, sparsity)
    else
        v.tmp_n .= abs.(v.β_prev .- v.β)
        return false, (maximum(v.tmp_n), sparsity)
    end
end

function cox_one_iter!(X::MPIArray, u::COXUpdate, v::COXVariables)
    copyto!(v.β_prev, v.β)
    update!(X, u, v)
end

function loop!(X::MPIArray, u, iterfun, evalfun, args...)
    converged = false
    t = 0
    while !converged && t < u.maxiter
        t += 1
        iterfun(X, u, args...)
        if t % u.step == 0
            converged, monitor = evalfun(X, u, args...)
            if DistStat.Rank() == 0
                println(t, ' ', monitor)
            end
        end
    end
end

function cox!(X::MPIArray, u::COXUpdate, v::COXVariables)
    loop!(X, u, cox_one_iter!, get_objective!, v)
end

include("cmdline.jl")
opts = parse_commandline_cox()
if DistStat.Rank() == 0
    println("world size: ", DistStat.Size())
    println(opts)
end

m = opts["rows"]
n = opts["cols"]
iter = opts["iter"]
interval = opts["step"]
T = Float64
A = Array
if opts["gpu"]
    using CUDA
    A = CuArray
end
if opts["Float32"]
    T = Float32
end
init_opt = opts["init_from_master"]
seed = opts["seed"]
eval_obj = opts["eval_obj"]
censor_rate = opts["censor_rate"]
lambda = opts["lambda"]

X = MPIMatrix{T, A}(undef, m, n)
rand!(X; common_init=init_opt, seed=seed)
δ = convert(A{T}, rand(m) .> censor_rate)
DistStat.Bcast!(δ) # synchronize the choice for δ.

uquick = COXUpdate(;maxiter=2, step=1, verbose=true)
u = COXUpdate(;maxiter=iter, step=interval, verbose=true)
# for simulation run, we just assume that the data are in reversed order of survivial time
t = convert(A{T}, collect(reverse(1:size(X,1))))
v = COXVariables(X, δ, lambda, t; eval_obj=eval_obj)
cox!(X, uquick, v)
reset!(v; seed=seed)
if DistStat.Rank() == 0
    @time cox!(X, u, v)
else
    cox!(X, u, v)
end


text:
export PrimesSet

################################################################################
#
#  Missing isprime functionality
#
################################################################################

# Fallback
function isprime(x::Integer)
  return isprime(fmpz(x))
end

function next_prime(x::BigInt, proved::Bool = true)
  return BigInt(next_prime(fmpz(x), proved))
end

function next_prime(x::T, proved::Bool = true) where {T <: Integer}
  return T(next_prime(BigInt(x), proved))
end

################################################################################
#
#  Primes iterator
#
################################################################################

# TODO (Tommy):
# At the moment proof doesn't do anything. The only reason is that
# next_prime() and isprime() do not support it generically.
# Once they do, it is easy to fix.
struct PrimesSet{T}
  from::T
  to::T
  nocond::Bool  # If set, there is no condition on the prime.
  proof::Bool   # Does not do anything at the moment
  mod::T        # If set (i.e. > 1), only primes p % mod == a are returned
  a::T
  sv::UInt

  function PrimesSet{T}(f::T, t::T) where {T}
    z = new{T}(f, t, true, false, T(1), T(0), UInt(1))
    return z
  end

  function PrimesSet{T}(f::T, t::T, mod::T, val::T) where {T}
    sv = UInt(1)
    p = UInt(2)
    while sv < 2^30 && p < f
      if !iszero(mod % p)
        sv *= p
      end
      p = next_prime(p, false)
    end
    mod = abs(mod)
    val = Base.mod(val, mod)
    if mod == 1 && val == 0
      return new{T}(f, t, false, false, mod, val, sv)
    end
    if mod < 2
      error("modulus needs to be > 1")
    end
    if val < 1 || val >= mod
      error("wrong congruence condition: value needs to be > 0 and < modulus")
    end
    if !isone(gcd(mod, val))
      error("Modulus and value need to be coprime.")
    end
    r = new{T}(f, t, false, false, mod, val, sv)
    return r
  end
end

@doc Markdown.doc"""
    PrimesSet(f::Integer, t::Integer) -> PrimesSet
    PrimesSet(f::fmpz, t::fmpz) -> PrimesSet

Returns an iterable object $S$ representing the prime numbers $p$
for $f \le p \le t$. If $t=-1$, then the upper bound is infinite.
"""
function PrimesSet(f::T, t::T) where T
  return PrimesSet{T}(f, t)
end

@doc Markdown.doc"""
    PrimesSet(f::Integer, t::Integer, mod::Integer, val::Integer)
    PrimesSet(f::fmpz, t::fmpz, mod::fmpz, val::fmpz)

Returns an iterable object $S$ representing the prime numbers $p$
for $f \le p \le t$ and $p\equiv val \bmod mod$ (primes in arithmetic
progression).
 If $t=-1$, then the upper bound is infinite.
"""
function PrimesSet(f::T, t::T, mod::T, val::T) where {T}
  return PrimesSet{T}(f, t, mod, val)
end

function Base.iterate(A::PrimesSet{T}) where {T <: Integer}
  if A.to != -1 && A.from > A.to
    return nothing
  end

  if A.nocond
    if !isprime(fmpz(A.from))
      p = next_prime(A.from)
    else
      p = A.from
    end

    if A.to != -1 && p > A.to
      return nothing
    else
      if p % A.mod != A.a
        return nothing
      else
        @show return p, p
      end
    end
  end

  curr = A.from
  c = curr % A.mod
  if A.mod > 1 && c != A.a
    curr += (- c + A.a) % A.mod
  end
  c_U = curr % A.sv
  c_M = A.mod % A.sv
  UIntone = one(UInt)
  i = zero(UInt)
  while gcd(c_U + i * c_M, A.sv) != UInt(1) || !isprime(fmpz(curr))
    curr += A.mod
    i += UIntone
    if A.to != -1 && curr > A.to
      break
    end
  end

  if A.to != -1 && curr > A.to
    return nothing
  else
    return curr, curr
  end
end

function Base.iterate(A::PrimesSet{fmpz})
  if A.to != -1 && A.from > A.to
    return nothing
  end

  if A.nocond
    if !isprime(A.from)
      p = next_prime(A.from)
    else
      p = A.from
    end

    if A.to != -1 && p > A.to
      return nothing
    else
      if p % A.mod != A.a
        return nothing
      else
        return p, p
      end
    end
  end

  curr = A.from
  c = curr % A.mod
  if A.mod > 1 && c != A.a
    curr += (-c + A.a) % A.mod
  end

  c_U = curr % A.sv
  c_M = A.mod % A.sv
  UIntone = one(UInt)
  i = zero(UInt)
  while !isone(gcd(c_U + i * c_M, A.sv)) || !isprime(fmpz(curr))
    curr += A.mod
    i += UIntone
    if A.to != -1 && curr > A.to
      break
    end
  end

  if A.to != -1 && curr > A.to
    return nothing
  else
    return curr, curr
  end
end

function Base.iterate(A::PrimesSet{T}, p) where T<: Union{Integer, fmpz}
  if A.nocond
    if p == 2
      nextp = T(3)
    else
      if T == Int
        nextp = next_prime(p, A.proof)
      else
        nextp = next_prime(p)
      end
    end

    if A.to != -1 && nextp > A.to
      return nothing
    else
      return nextp, nextp
    end
  end

  if A.mod > 1
    m = A.mod
  else
    if p==2
      nextp = T(3)
    end
    m = T(2)
  end
  nextp = p+m
  i = zero(UInt)
  UIntone = one(UInt)
  c_U = nextp % A.sv
  c_M = m % A.sv
  while !isone(gcd(c_U + i * c_M, A.sv)) || !isprime(fmpz(nextp))
    nextp += m
    i += UIntone
    if A.to != -1 && nextp > A.to
      break
    end
  end

  if A.to != -1 && nextp > A.to
    return nothing
  else
    return nextp, nextp
  end
end

# Iteration interface
#function Base.start(A::PrimesSet{T}) where T <: Integer
#  if A.nocond
#    if !isprime(fmpz(A.from))
#      p = next_prime(A.from)
#    else
#      p = A.from
#    end
#    return p
#  end
#
#  curr = A.from
#  c = curr % A.mod
#  if A.mod > 1 && c != A.a
#    curr += (- c + A.a) % A.mod
#  end
#  c_U = curr % A.sv
#  c_M = A.mod % A.sv
#  UIntone = one(UInt)
#  i = zero(UInt)
#  while gcd(c_U + i * c_M, A.sv) != UInt(1) || !isprime(fmpz(curr))
#    curr += A.mod
#    i += UIntone
#  end
#  return curr
#end

#function Base.start(A::PrimesSet{fmpz})
#  if A.nocond
#    if !isprime(A.from)
#      p = next_prime(A.from)
#    else
#      p = A.from
#    end
#    return p
#  end
#
#  curr = A.from
#  c = curr % A.mod
#  if A.mod > 1 && c != A.a
#    curr += (-c + A.a) % A.mod
#  end
#
#  c_U = curr % A.sv
#  c_M = A.mod % A.sv
#  UIntone = one(UInt)
#  i = zero(UInt)
#  while !isone(gcd(c_U + i * c_M, A.sv)) || !isprime(fmpz(curr))
#    curr += A.mod
#    i += UIntone
#  end
#  return curr
#end
#
#function Base.next(A::PrimesSet{T}, st::T) where T<: Union{Integer, fmpz}
#  p = st
#  if A.nocond
#    if p == 2
#      st = T(3)
#      return p, st
#    else
#      if T == Int
#        st = next_prime(p, A.proof)
#      else
#        st = next_prime(p)
#      end
#      return p, st
#    end
#  end
#
#  if A.mod > 1
#    m = A.mod
#  else
#    if p==2
#      st = T(3)
#      return p, st
#    end
#    m = T(2)
#  end
#  st = p+m
#  i = zero(UInt)
#  UIntone = one(UInt)
#  c_U = st % A.sv
#  c_M = m % A.sv
#  while !isone(gcd(c_U + i * c_M, A.sv)) || !isprime(fmpz(st))
#    st += m
#    i += UIntone
#  end
#
#  return p, st
#end
#
#function Base.done(A::PrimesSet{T}, st::T) where T <: Union{Integer, fmpz}
#  return A.to != -1 && st > A.to
#end

Base.eltype(::PrimesSet{T}) where {T <: Union{Integer, fmpz}} = T

Base.length(A::PrimesSet) = length(collect(A))

Base.IteratorSize(::Type{PrimesSet{T}}) where {T} = Base.SizeUnknown()

#Base.iteratorsize(::Type{PrimesSet{T}}) where {T} = Base.SizeUnknown()


text:
using ModuleInfo
using Test
using ReTest

ModuleInfo.runtests()

@testset "ModuleInfo.jl" begin
    # Write your tests here.
end


text:
"""
    function calc_rate(su, rel_reserve::Tuple, turnover::Tuple,
                       j_E_mai, y_E_Ea, y_E_Eb, y_V_E, κsoma, tstep)

Calculate growth rate using a numeric root-finder, also determining
wether the organ is alive or dead.

Returns a `Tuple` holding the rate and a Bool for alive/dead status
"""
function calc_rate end
function calc_rate(su, rel_reserve::Tuple, turnover::Tuple,
                   j_E_mai, y_E_Ea, y_E_Eb, y_V_E, κsoma, tstep)
    # Find the type of `rate` so that diff and units work with atol etc
    one_r = oneunit(eltype(turnover))
    # Get rate with a zero finder
    r, state = let turnover=turnover, j_E_mai=j_E_mai, y_E_Ea=y_E_Ea, y_E_Eb=y_E_Eb, y_V_E=y_V_E, κsoma=κsoma
        atol = one_r * 1e-10; maxiters = 200
        findzero(Secant(-2one_r, 1one_r), atol, maxiters) do r
            rate_formula(r, su, rel_reserve, turnover, j_E_mai, y_E_Ea, y_E_Eb, y_V_E, κsoma)
        end
    end

    if state == DEAD
        @warn "Root for rate not found at t = $tstep"
        return zero(r), DEAD
    elseif r < zero(r)
        @warn "Rate is less than zero at t = $tstep"
        return zero(r), DEAD
    else
        return r, ALIVE
    end
end

"""
    rate_formula(r, ureserve::Tuple, turnover::Tuple, j_E_mai, y_V_E, κsoma)

Rate formulas for E, CN or CNE reserves
"""
function rate_formula end
rate_formula(r, su, rel_reserve::NTuple{1}, turnover::NTuple{1},
             j_E_mai, y_E_Ea, y_E_Eb, y_V_E, κsoma) = begin
    (j_E,) = rel_reserve .* (turnover .- r)
    y_V_E * (κsoma * j_E - j_E_mai) - r
end
rate_formula(r, su, rel_reserve::NTuple{2}, turnover::NTuple{2},
             j_E_mai, y_E_Ea, y_E_Eb, y_V_E, κsoma) = begin
    j_Ea, j_Eb = rel_reserve .* (turnover .- r)
    j_E = synthesizing_unit(su, j_Ea * y_E_Ea, j_Eb * y_E_Eb)
    y_V_E * (κsoma * j_E - j_E_mai) - r
end
rate_formula(r, su, rel_reserve::NTuple{3}, turnover::NTuple{3},
             j_E_mai, y_E_Ea, y_E_Eb, y_V_E, κsoma) = begin
    j_Ea, j_Eb, j_E = rel_reserve .* (turnover .- r)
    j_E += synthesizing_unit(su, j_Ea * y_E_Ea, j_Eb * y_E_Eb)
    y_V_E * (κsoma * j_E - j_E_mai) - r
end


text:
module BaseExtensions

export anynonzero
export allzero
export zero!
export realtype
export reinterpretreal
export reinterpretcomplex
export squeeze
export flipdims
export fillcopy
export fillcopy!
export filldeepcopy
export filldeepcopy!

"""
	anynonzero(a)

Check if `a` contains any non-zero element.

# Examples
```jldoctest
julia> a = [1, 2, 0];

julia> anynonzero(a)
true

julia> a = zeros(2);

julia> anynonzero(a)
false
```
"""
anynonzero(a) = any(!iszero, a)

"""
	allzero(a)

Check if all elements of `a` are zero.

# Examples
```jldoctest
julia> a = [1, 2, 0];

julia> allzero(a)
false

julia> a = zeros(2);

julia> allzero(a)
true
```
"""
allzero(a) = all(iszero, a)

"""
	zero!(a::AbstractArray)

Set all elements of an `AbstractArray` to zero.

# Examples

```jldoctest
julia> a = [1, 2];

julia> zero!(a);

julia> a
2-element $(Array{Int,1}):
 0
 0
```
"""
zero!(a::AbstractArray) = fill!(a, zero(eltype(a)))

"""
	realtype(::Complex)
	realtype(::Type{<:Complex})

Return the type of the real and imaginary parts of a `Complex` number.

# Examples

```jldoctest
julia> realtype(ComplexF64(2,2))
Float64

julia> realtype(ComplexF32(2,2))
Float32

julia> realtype(ComplexF32)
Float32
```
"""
realtype(::Type{Complex{T}}) where {T} = T
realtype(c::Complex) = realtype(typeof(c))

"""
	reinterpretreal(arr::AbstractArray{<:Complex})

Reinterpret a `Complex` array as a `Real` one containing the real and imaginary parts.

# Examples
```jldoctest
julia> a = [1 + 3im, 2+ 4im]
2-element $(Array{Complex{Int},1}):
 1 + 3im
 2 + 4im

julia> reinterpretreal(a)
4-element reinterpret($Int, ::$(Array{Complex{Int},1})):
 1
 3
 2
 4
```

See also: [`reinterpretcomplex`](@ref)
"""
function reinterpretreal(arr::AbstractArray{<:Complex})
	reinterpret(realtype(eltype(arr)), arr)
end

"""
	reinterpretcomplex(arr::AbstractArray{<:Real})

Reinterpret a `Real` array as a `Complex` one assuming that the real and imaginary parts
are laid out consecutively.

# Examples
```jldoctest
julia> a = [1, 2, 3, 4];

julia> reinterpretcomplex(a)
2-element reinterpret(Complex{$Int}, ::$(Array{Int,1})):
 1 + 2im
 3 + 4im
```

See also: [`reinterpretreal`](@ref)
"""
function reinterpretcomplex(arr::AbstractArray{<:Real})
	reinterpret(Complex{eltype(arr)}, arr)
end

"""
	squeeze(A)

Remove dimensions from `A` that have `size(A, dim) == 1`.
Returns a view that shares the underlying data with `A`.

# Examples
```jldoctest
julia> A = reshape(collect(1:4), 1, 4)
1×4 $(Array{Int,2}):
 1  2  3  4

julia> B = squeeze(A)
4-element $(Array{Int,1}):
 1
 2
 3
 4
```
"""
function squeeze(A)
	dims = Tuple(findall(==(1), size(A)))
	dropdims(A, dims = dims)
end

"""
	flipdims(A)

Permute `A` so that the sequence of axes is reversed. For 2D arrays this is an eager transpose.
The operation is not recursive.

# Examples
```jldoctest
julia> A = reshape(collect(1:8), 2, 4)
2×4 $(Array{Int,2}):
 1  3  5  7
 2  4  6  8

julia> B = flipdims(A)
4×2 $(Array{Int,2}):
 1  2
 3  4
 5  6
 7  8

julia> A = reshape(collect(1:8), 2, 2, 2)
2×2×2 $(Array{Int,3}):
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 5  7
 6  8

julia> flipdims(A)
2×2×2 $(Array{Int,3}):
[:, :, 1] =
 1  3
 5  7

[:, :, 2] =
 2  4
 6  8
```
"""
function flipdims(A)
	permutedims(A, reverse(collect(1:ndims(A))))
end

"""
	fillcopy(x, dims...)

Create an array filled with shallow copies of `x`, where the dimensions or indices of the array is
specfied through the keyword `dims`. As in `copy`, only the outer container of `x` is copied.

Unlike `fill`, returns an array where each element is unique.

# Examples
```jldoctest
julia> A = fillcopy(ones(2), 2, 2)
2×2 $(Array{Array{Float64,1},2}):
 [1.0, 1.0]  [1.0, 1.0]
 [1.0, 1.0]  [1.0, 1.0]

julia> A[1,1] === A[2,2]
false

julia> A[1,1] .= 2;

julia> A
2×2 $(Array{Array{Float64,1},2}):
 [2.0, 2.0]  [1.0, 1.0]
 [1.0, 1.0]  [1.0, 1.0]

julia> B = fillcopy(A, 2)
2-element $(Array{Array{Array{Float64,1},2},1}):
 [[2.0, 2.0] [1.0, 1.0]; [1.0, 1.0] [1.0, 1.0]]
 [[2.0, 2.0] [1.0, 1.0]; [1.0, 1.0] [1.0, 1.0]]

julia> B[2] === B[1]
false

julia> B[2][1] === B[1][1]
true
```

See also: [`filldeepcopy`](@ref)
"""
function fillcopy(x, inds...)
	A = similar(Array{typeof(x)}, inds)
	fillcopy!(A, x)
	return A
end

"""
	fillcopy!(A, x)

Fill `A` with shallow copies of `x`. As in `copy`, only the outer container of `x` is copied.

# Examples
```jldoctest
julia> A = Vector{Vector{Vector{Int}}}(undef, 3);

julia> fillcopy!(A, [[1,2], [1, 2,3]])
3-element $(Array{Array{Array{Int,1},1},1}):
 [[1, 2], [1, 2, 3]]
 [[1, 2], [1, 2, 3]]
 [[1, 2], [1, 2, 3]]

julia> A[1][1] === A[2][1]
true

julia> A[1][1] .= 4:5;

julia> A
3-element $(Array{Array{Array{Int,1},1},1}):
 [[4, 5], [1, 2, 3]]
 [[4, 5], [1, 2, 3]]
 [[4, 5], [1, 2, 3]]
```

See also: [`filldeepcopy!`](@ref)
"""
function fillcopy!(A, x)
	@inbounds for i in eachindex(A)
		A[i] = copy(x)
	end
	return A
end

"""
	filldeepcopy(x, dims...)

Create an array filled with deep copies of `x`, where the dimensions or indices of the array is
specfied through the keyword `dims`. As in `deepcopy`, the copy operation is recursive, and each element
in the returned array will be completely unique in identity, albeit equal to `x` in value.

# Examples
```jldoctest
julia> A = filldeepcopy(ones(2), 2, 2)
2×2 $(Array{Array{Float64,1},2}):
 [1.0, 1.0]  [1.0, 1.0]
 [1.0, 1.0]  [1.0, 1.0]

julia> A[1,1] === A[2,2]
false

julia> A[1,1] .= 2;

julia> A
2×2 $(Array{Array{Float64,1},2}):
 [2.0, 2.0]  [1.0, 1.0]
 [1.0, 1.0]  [1.0, 1.0]

julia> B = filldeepcopy(A, 2)
2-element $(Array{Array{Array{Float64,1},2},1}):
 [[2.0, 2.0] [1.0, 1.0]; [1.0, 1.0] [1.0, 1.0]]
 [[2.0, 2.0] [1.0, 1.0]; [1.0, 1.0] [1.0, 1.0]]

julia> B[2] === B[1]
false

julia> B[2][1] === B[1][1]
false
```

See also: [`fillcopy`](@ref)
"""
function filldeepcopy(x, inds...)
	A = similar(Array{typeof(x)}, inds)
	filldeepcopy!(A, x)
	return A
end

"""
	filldeepcopy!(A, x)

Fill `A` with deep copies of `x`. The operation is recursive and each element in `A` will be
completely unique in identity, albeit equal to `x` in value.

# Examples
```jldoctest
julia> A = Vector{Vector{Vector{Int}}}(undef, 3);

julia> filldeepcopy!(A, [[1,2], [1, 2,3]])
3-element $(Array{Array{Array{Int,1},1},1}):
 [[1, 2], [1, 2, 3]]
 [[1, 2], [1, 2, 3]]
 [[1, 2], [1, 2, 3]]

julia> A[1][1] === A[2][1]
false

julia> A[1][1] .= 4:5;

julia> A
3-element $(Array{Array{Array{Int,1},1},1}):
 [[4, 5], [1, 2, 3]]
 [[1, 2], [1, 2, 3]]
 [[1, 2], [1, 2, 3]]
```

See also: [`fillcopy!`](@ref)
"""
function filldeepcopy!(A, x)
	@inbounds for i in eachindex(A)
		A[i] = deepcopy(x)
	end
	return A
end

end


text:
function _move(nrows::Int, ncols::Int, cell::Cell, x::Int, y::Int)
    col = (cell.col + x - 1 + ncols) % ncols + 1
    row = (cell.row + y - 1 + nrows) % nrows + 1
    return row, col
end

function _move(grid::Grid, cell::Cell, x, y)
    row, col = _move(grid.nrows, grid.ncols, cell, x, y)
    return grid[row, col]
end

north(nrows, ncols, cell::Cell) = _move(nrows, ncols, cell, 0, -1)
south(nrows, ncols, cell::Cell) = _move(nrows, ncols, cell, 0, 1)
east(nrows, ncols, cell::Cell{Rect}) = _move(nrows, ncols, cell, 1, 0)
west(nrows, ncols, cell::Cell{Rect}) = _move(nrows, ncols, cell, -1, 0)
east(nrows, ncols, cell::Cell{Tria}) = _move(nrows, ncols, cell, 1, 0)
west(nrows, ncols, cell::Cell{Tria}) = _move(nrows, ncols, cell, -1, 0)

north_diagonal_y(cell::Cell{Hexa}) = iseven(cell.col) ? -1 : 0
south_diagonal_y(cell::Cell{Hexa}) = iseven(cell.col) ? 0 : 1

function northwest(nrows, ncols, cell::Cell{Hexa})
    return _move(nrows, ncols, cell, -1, north_diagonal_y(cell))
end
function northeast(nrows, ncols, cell::Cell{Hexa})
    return _move(nrows, ncols, cell, 1, north_diagonal_y(cell))
end
function southwest(nrows, ncols, cell::Cell{Hexa})
    return _move(nrows, ncols, cell, -1, south_diagonal_y(cell))
end
function southeast(nrows, ncols, cell::Cell{Hexa})
    return _move(nrows, ncols, cell, 1, south_diagonal_y(cell))
end

north(grid::Grid, cell::Cell) = _move(grid, cell, 0, -1)
south(grid::Grid, cell::Cell) = _move(grid, cell, 0, 1)
east(grid::Grid{Rect}, cell::Cell{Rect}) = _move(grid, cell, 1, 0)
west(grid::Grid{Rect}, cell::Cell{Rect}) = _move(grid, cell, -1, 0)
east(grid::Grid{Tria}, cell::Cell{Tria}) = _move(grid, cell, 1, 0)
west(grid::Grid{Tria}, cell::Cell{Tria}) = _move(grid, cell, -1, 0)

function northwest(grid::Grid{Hexa}, cell::Cell{Hexa})
    return _move(grid, cell, -1, north_diagonal_y(cell))
end
function northeast(grid::Grid{Hexa}, cell::Cell{Hexa})
    return _move(grid, cell, 1, north_diagonal_y(cell))
end
function southwest(grid::Grid{Hexa}, cell::Cell{Hexa})
    return _move(grid, cell, -1, south_diagonal_y(cell))
end
function southeast(grid::Grid{Hexa}, cell::Cell{Hexa})
    return _move(grid, cell, 1, south_diagonal_y(cell))
end

north(grid::Grid, row, col) = north(grid, grid[row, col])
south(grid::Grid, row, col) = south(grid, grid[row, col])
east(grid::Grid{Rect}, row, col) = east(grid, grid[row, col])
west(grid::Grid{Rect}, row, col) = west(grid, grid[row, col])
east(grid::Grid{Tria}, row, col) = east(grid, grid[row, col])
west(grid::Grid{Tria}, row, col) = west(grid, grid[row, col])

northwest(grid::Grid{Hexa}, row, col) = northwest(grid, grid[row, col])
northeast(grid::Grid{Hexa}, row, col) = northeast(grid, grid[row, col])
southwest(grid::Grid{Hexa}, row, col) = southwest(grid, grid[row, col])
southeast(grid::Grid{Hexa}, row, col) = southeast(grid, grid[row, col])


text:
module JSONPointer

using OrderedCollections

include("pointer.jl")

export @j_str

end # module


text:
export MPI_rank, MPI_size

MPI_rank = 0
MPI_size = 1


text:
function cal_flux_residual_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, numPoints)
	idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
	thread_idx = threadIdx().x
	block_dim = blockDim().x

	shared = @cuStaticSharedMem(Float64, 1024)
	flux_shared = @cuStaticSharedMem(Float64, 512)
	if idx > 0 && idx <= numPoints

		flux_shared[thread_idx], flux_shared[thread_idx + block_dim], flux_shared[thread_idx + block_dim * 2],
			flux_shared[thread_idx + block_dim * 3] = 0.0, 0.0, 0.0, 0.0

		if gpuGlobalDataFixedPoint[idx].flag_1 == gpuConfigData[17]
			wall_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)
		elseif gpuGlobalDataFixedPoint[idx].flag_1 == gpuConfigData[18]
			interior_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)
		elseif gpuGlobalDataFixedPoint[idx].flag_1 == gpuConfigData[19]
			outer_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)
		end
	else
		return nothing
	end

	gpuGlobalDataRest[5, idx] = flux_shared[thread_idx]
	gpuGlobalDataRest[6, idx] = flux_shared[thread_idx + block_dim]
	gpuGlobalDataRest[7, idx] = flux_shared[thread_idx + block_dim * 2]
	gpuGlobalDataRest[8, idx] = flux_shared[thread_idx + block_dim * 3]
	sync_threads()
	return nothing
end

@inline function wall_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)
	wall_dGx_pos_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)

	wall_dGx_neg_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)

	wall_dGy_neg_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)
	return nothing
end

@inline function interior_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)

	interior_dGx_pos_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)

	interior_dGx_neg_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)

	interior_dGy_pos_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)

	interior_dGy_neg_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)
	# interior_dG_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)
	return nothing
end

@inline function outer_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)
	outer_dGx_pos_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)
	outer_dGx_neg_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)
	outer_dGy_pos_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)
	return nothing
end

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
#=
    Model Predictive Control (feedback) approach for a mechanical oscillator with two masses, springs and attenuators

    Inspired by: https://github.com/rdeits/DynamicWalking2018.jl

    Author: Stephan Scholz
    Year: 2021
=#

using JuMP, Ipopt

my_optimizer = with_optimizer(Ipopt.Optimizer, print_level = 0)

const d₁ = 0.5; # Damping
const d₂ = 0.5;
const c₁ = 0.1; # Spring constant
const c₂ = 0.1;
const m₁ = 1.0; # Mass
const m₂ = 2.0;

# System matrix
const A = [0 0 1 0; 0 0 0 1; (-1/m₁)*(c₁ + c₂) c₂/m₁ (-1/m₁)*(d₁ + d₂) d₂/m₁; c₂/m₂ -c₂/m₂ d₂/m₂ -d₂/m₂]

const B = [0 0; 0 0; 1/m₁ 0; 0 1/m₂]  # Input matrix
const C = [1 0 0 0; 0 1 0 0]          # Output matrix

mech_osc(x,u) = A * x + B * u   # System dynamics

function mpc_emo(x_init)

    model = Model(my_optimizer)

    ΔT = 0.1            #  Sampling time
    n_horizon = 20;     # MPC horizon

    Nx = size(A)[1]
    Nu = size(B)[2]

    # Input constraints
    umin = -100.0
    umax = 100.0

    ref = [2.0, 4.0, 0.0, 0.0] # Reference

    @variables model begin
        x[1:Nx, 1:n_horizon]    # States
        e[1:Nx, 1:n_horizon]    # Errors = Reference - States
        umin <= u[1:Nu, 1:n_horizon] <= umax  # Input constraints
    end

    @constraint(model, x[1:Nx, 1] .== x_init )                          # Initial values
    @constraint(model, [k=1:n_horizon], e[:, k] .== ref[:] - x[:, k])   # Error = Reference - States
    @constraint(model, [k=1:n_horizon-1], x[:, k+1] .== x[:, k] + ΔT * mech_osc(x[:,k], u[:,k])) # System dynamics

    J = @NLexpression(model, 100*sum(e[i,end]^2 for i in 1:Nx) + sum(u[j,end]^2 for j in 1:Nu) )

    @NLobjective( model, Min, J)

    optimize!(model)
    return JuMP.value.(x), JuMP.value.(u), JuMP.value.(e)
end

xpos = rand(4)          # Recent states
x_hist = zeros(4,0)     # History of states

for i in 1:50

    # Run the MPC control optimization
    x_plan, u_plan, e_plan = mpc_emo(xpos)

    # Save states
    x_hist = hcat(x_hist, x_plan[:,1])

    # Apply the planned input signals and simulate one step in time
    ΔT = 0.1
    xpos = xpos + ΔT * mech_osc(xpos, u_plan[:,1])
end

using Plots
plot(x_hist', labels=["x1" "x2" "x3" "x4"], title="Model Predictive Control of a mechanical oscillator")


text:
# This example requires Cloudy to be installed.
#using Pkg; Pkg.add(PackageSpec(name="Cloudy", version="0.1.0"))
using Cloudy
const PDistributions = Cloudy.ParticleDistributions

# Import modules
using Distributions  # probability distributions and associated functions
using StatsBase
using LinearAlgebra
using StatsPlots
using Plots
using JLD2
using Random

# Import Calibrate-Emulate-Sample modules
using EnsembleKalmanProcesses.EnsembleKalmanProcessModule
using EnsembleKalmanProcesses.Observations
using EnsembleKalmanProcesses.ParameterDistributionStorage

# Import the module that runs Cloudy
include("GModel.jl")
using .GModel

################################################################################
#                                                                              #
#                      Cloudy Calibrate-Emulate-Sample Example                 #
#                                                                              #
#                                                                              #
#     This example uses Cloudy, a microphysics model that simulates the        #
#     coalescence of cloud droplets into bigger drops, to demonstrate how      #
#     the full Calibrate-Emulate-Sample pipeline can be used for Bayesian      #
#     learning and uncertainty quantification of parameters, given some        #
#     observations.                                                            #
#                                                                              #
#     Specifically, this examples shows how to learn parameters of the         #
#     initial cloud droplet mass distribution, given observations of some      #
#     moments of that mass distribution at a later time, after some of the     #
#     droplets have collided and become bigger drops.                          #
#                                                                              #
#     In this example, Cloudy is used in a "perfect model" (aka "known         #
#     truth") setting, which means that the "observations" are generated by    #
#     Cloudy itself, by running it with the true parameter values. In more     #
#     realistic applications, the observations will come from some external    #
#     measurement system.                                                      #
#                                                                              #
#     The purpose is to show how to do parameter learning using                #
#     Calibrate-Emulate-Sample in a simple (and highly artificial) setting.    #
#                                                                              #
#     For more information on Cloudy, see                                      #
#              https://github.com/CliMA/Cloudy.jl.git                          #
#                                                                              #
################################################################################

rng_seed = 41
Random.seed!(rng_seed)

homedir = pwd()
println(homedir)
figure_save_directory = homedir*"/output/"
data_save_directory = homedir*"/output/"
if ~isdir(figure_save_directory)
    mkdir(figure_save_directory)
end
if ~isdir(data_save_directory)
    mkdir(data_save_directory)
end

###
###  Define the (true) parameters and their priors
###

# Define the parameters that we want to learn
# We assume that the true particle mass distribution is a Gamma distribution
# with parameters N0_true, θ_true, k_true
param_names = ["N0", "u", "k"]
n_param = length(param_names)

N0_true = 300.0  # number of particles (scaling factor for Gamma distribution)
θ_true = 1.5597  # scale parameter of Gamma distribution
k_true = 0.0817  # shape parameter of Gamma distribution
params_true = [N0_true, θ_true, k_true]
# Note that dist_true is a Cloudy distribution, not a Distributions.jl
# distribution
dist_true = PDistributions.Gamma(N0_true, θ_true, k_true)

###
###  Define priors for the parameters we want to learn
###

# Define constraints
lbound_N0 = 0.4 * N0_true
lbound_θ = 1.0e-1
lbound_k = 1.0e-4
c1 = bounded_below(lbound_N0)
c2 = bounded_below(lbound_θ)
c3 = bounded_below(lbound_k)
constraints = [[c1], [c2], [c3]]

# We choose to use normal distributions to represent the prior distributions of
# the parameters in the transformed (unconstrained) space. i.e log coordinates
d1 = Parameterized(Normal(4.5, 1.0)) #truth is 5.19
d2 = Parameterized(Normal(0.0, 2.0)) #truth is 0.378
d3 = Parameterized(Normal(-1.0, 1.0))#truth is -2.51
distributions = [d1, d2, d3]

prior_mean = [4.5; 0.0; -1.0]
prior_cov = Array(Diagonal( ones(Float64, length(prior_mean) ) ) )

param_names = ["N0", "u", "k"]

priors = ParameterDistribution(distributions, constraints, param_names)

###
###  Define the data from which we want to learn the parameters
###

data_names = ["M0", "M1", "M2"]
moments = [0.0, 1.0, 2.0]
n_moments = length(moments)

###
###  Model settings
###

# Collision-coalescence kernel to be used in Cloudy
coalescence_coeff = 1/3.14/4/100
kernel_func = x -> coalescence_coeff
kernel = Cloudy.KernelTensors.CoalescenceTensor(kernel_func, 0, 100.0)

# Time period over which to run Cloudy
tspan = (0., 1.0)

###
###  Generate (artificial) truth samples
###  Note: The observables y are related to the parameters u by:
###        y = G(x1, x2) + η
###

g_settings_true = GModel.GSettings(kernel, dist_true, moments, tspan)
gt = GModel.run_G(params_true, g_settings_true, PDistributions.update_params,
                  PDistributions.moment, Cloudy.Sources.get_int_coalescence)

n_samples = 100
yt = zeros(length(gt),n_samples)
# In a perfect model setting, the "observational noise" represent the internal
# model variability. Since Cloudy is a purely deterministic model, there is no
# straightforward way of coming up with a covariance structure for this internal
# model variability. We decide to use a diagonal covariance, with entries
# (variances) largely proportional to their corresponding data values, gt.
Γy = convert(Array, Diagonal([100.0, 5.0, 30.0]))
μ = zeros(length(gt))

# Add noise
for i in 1:n_samples
    yt[:, i] = gt .+ rand(MvNormal(μ, Γy))
end

truth = Observations.Obs(yt, Γy, data_names)
truth_sample = truth.mean

###
###  Calibrate: Unscented Kalman Inversion
###

N_iter = 20 # number of iterations
# need to choose regularization factor α ∈ (0,1],
# when you have enough observation data α=1: no regularization
α_reg =  1.0
# update_freq 1 : approximate posterior covariance matrix with an uninformative prior
#             0 : weighted average between posterior covariance matrix with an uninformative prior and prior
update_freq = 1

process = Unscented(prior_mean, prior_cov, length(truth_sample), α_reg, update_freq)
ukiobj = EnsembleKalmanProcessModule.EnsembleKalmanProcess(truth_sample, truth.obs_noise_cov, process)

# Initialize a ParticleDistribution with dummy parameters. The parameters
# will then be set in run_G_ensemble
dummy = 1.0
dist_type = PDistributions.Gamma(dummy, dummy, dummy)
g_settings = GModel.GSettings(kernel, dist_type, moments, tspan)

err = zeros(N_iter)
for i in 1:N_iter

    params_i = mapslices(x -> transform_unconstrained_to_constrained(priors, x),
                         get_u_final(ukiobj); dims=1)

    # # prediction step
    # u_p = EnsembleKalmanProcessModule.update_ensemble_prediction!(ukiobj)

    # # define black box parameter to observation map,
    # # with certain parameter transformation related to imposing some constraints
    # # i.e. θ -> e^θ  -> G(e^θ) = y
    # u_p_trans = mapslices(x -> transform_unconstrained_to_constrained(priors, x), u_p; dims=1)
    g_ens = GModel.run_G_ensemble(params_i, g_settings,
                                  PDistributions.update_params,
                                  PDistributions.moment,
                                  Cloudy.Sources.get_int_coalescence)
    # analysis step
    EnsembleKalmanProcessModule.update_ensemble!(ukiobj, g_ens)

    err[i] = get_error(ukiobj)[end] #mean((params_true - mean(params_i,dims=2)).^2)
    println("Iteration: "*string(i)*", Error: "*string(err[i])*" norm(Cov): "*string(norm(ukiobj.process.uu_cov[i])))
end

# UKI results: the mean is in ukiobj.process.u_mean
#              the covariance matrix is in ukiobj.process.uu_cov
transformed_params_true = transform_constrained_to_unconstrained(priors,
                                                                 params_true)

println("True parameters (transformed): ")
println(transformed_params_true)

println("\nUKI results:")
println(get_u_mean_final(ukiobj))

####
θ_mean_arr = hcat(ukiobj.process.u_mean...)
N_θ = length(ukiobj.process.u_mean[1])
θθ_std_arr = zeros(Float64, (N_θ, N_iter+1))
for i = 1:N_iter+1
    for j = 1:N_θ
        θθ_std_arr[j, i] = sqrt(ukiobj.process.uu_cov[i][j,j])
    end
end

ites = Array(LinRange(1, N_iter+1, N_iter+1))
plot(ites,grid=false, θ_mean_arr[1,:], yerror=3.0*θθ_std_arr[1,:],  label="u1")
plot!(ites, fill(transformed_params_true[1], N_iter+1), linestyle=:dash, linecolor=:grey,label=nothing)
plot!(ites,grid=false, θ_mean_arr[2,:], yerror=3.0*θθ_std_arr[2,:], label="u2", xaxis="Iterations")
plot!(ites, fill(transformed_params_true[2], N_iter+1), linestyle=:dash, linecolor=:grey,label=nothing)
plot!(ites,grid=false, θ_mean_arr[3,:], yerror=3.0*θθ_std_arr[3,:], label="u3", xaxis="Iterations")
plot!(ites, fill(transformed_params_true[3], N_iter+1), linestyle=:dash, linecolor=:grey,label=nothing)

#plots
gr(size=(1800,600))

for i in 1:N_iter
    θ_mean, θθ_cov = ukiobj.process.u_mean[i], ukiobj.process.uu_cov[i]
    θ1, θ2, fθ1θ2 = Gaussian_2d(θ_mean[1:2], θθ_cov[1:2,1:2], 100, 100)
    p1 = contour(θ1, θ2, fθ1θ2)
    plot!(p1,[transformed_params_true[1]], xaxis="u1", yaxis="u2", seriestype="vline",
        linestyle=:dash, linecolor=:red, label = false,
        title = "UKI iteration = " * string(i)
        )
    plot!(p1,[transformed_params_true[2]], seriestype="hline", linestyle=:dash, linecolor=:red, label = "optimum")

    θ2, θ3, fθ2θ3 = Gaussian_2d(θ_mean[2:3], θθ_cov[2:3,2:3], 100, 100)
    p2 = contour(θ2, θ3, fθ2θ3)
    plot!(p2,[transformed_params_true[2]], xaxis="u2", yaxis="u3", seriestype="vline",
        linestyle=:dash, linecolor=:red, label = false,
        title = "UKI iteration = " * string(i)
        )
    plot!(p2,[transformed_params_true[3]], seriestype="hline", linestyle=:dash, linecolor=:red, label = "optimum")

    θ3, θ1, fθ3θ1 = Gaussian_2d(θ_mean[3:-2:1], θθ_cov[3:-2:1,3:-2:1], 100, 100)
    p3 = contour(θ3, θ1, fθ3θ1)
    plot!(p3,[transformed_params_true[3]], xaxis="u3", yaxis="u1", seriestype="vline",
        linestyle=:dash, linecolor=:red, label = false,
        title = "UKI iteration = " * string(i)
        )
    plot!(p3,[transformed_params_true[1]], seriestype="hline", linestyle=:dash, linecolor=:red, label = "optimum")

    p=plot(p1, p2, p3, layout=(1,3))
    display(p)
    sleep(0.5)
end


text:
using CategoricalArrays
import Statistics

"""
    reorder!(fac::::CategoricalVector, v::Vector{<:Real}, f::Function=Statistics.mean)

Assign the levels of the ordered categorical vector `fac` in increasing order of `f(v)`
applied to what would be written in R as `split(v)`.
"""
function reorder!(fac::CategoricalVector, v::Vector{<:Real}, f::Function=Statistics.mean)
    if length(fac) ≠ length(v)
        throw(DimensionMismatch("lengths of fac and v must be equal"))
    end
    dict = Dict(l => eltype(v)[] for l in levels(fac))
    for (l, vv) in zip(fac, v)
        push!(dict[l], vv)
    end
    levels!(fac, collect(keys(sort(Dict(k => f(v) for (k,v) in dict), byvalue=true))))
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

module Utils

using StaticArrays
using LinearAlgebra: norm
using JuLIP: JVecF
using ACE: TransformedJacobi,  inv_transform,
            PolyTransform, PolyCutoff1s, SparseSHIP, SHIPBasis,
            TransformedPolys

import ACE
import ACE: rand_sphere, rand_radial
import Base: rand
import JuLIP: evaluate
import JuLIP.MLIPs: IPBasis

_get_ll(KL, νz) = getfield.(KL[νz.ν], :l)

# Some utilities to analyze the basis
# ------------------------------------

"""
quickly  generate a basis for testing and analysing
"""
function TestBasis(N, deg; rnn=1.0, rcut = 2.5, rin = 0.5, kwargs... )
   trans = PolyTransform(2, rnn)
   fcut = PolyCutoff1s(2, rin, rcut)
   spec = SparseSHIP(N, deg)
   return SHIPBasis(spec, trans, fcut; kwargs...)
end

"""
find all basis functions with prescribed ll tuple

TODO: generalize so we can also prescribe kk, or both.
"""

function findall_basis(basis; N = nothing, ll = nothing,
                              kwargs...)
   if N != nothing &&  ll != nothing
      error("findall_basis: only tell me N or ll but not both")
   end
   if N != nothing
      return findall_basis_N(basis, N; kwargs...)
   end
   if ll != nothing
      return findall_basis_ll(basis, ll; kwargs...)
   end
   error("findall_basis: I need either an N or an ll argument")
end

function findall_basis_N(basis, N; verbose=true)::Vector{Int}
   verbose && @warn("This code assumes there is only one species. (not checked!)")
   I1 = findall(b -> length(b[end]) == N, basis.bgrps[1])
   if isempty(I1); return Int[]; end
   verbose && @info("There are $(length(I1)) basis groups with N = $N:")
   Ib = mapreduce(i -> (basis.firstb[1][i]+1):basis.firstb[1][i+1],
                  vcat, I1)
   verbose && @info("  ... and $(length(Ib)) basis functions.")
   return Ib
end

function findall_basis_ll(basis, ll)
   ll = SVector(ll...)
   @warn("This code assumes there is only one species. (not checked!)")
   @info(" ll = $(ll)")
   @info("Get the purely rotation-invariant basis:")
   CA = ACE.Rotations.Rot3DCoeffs()
   Brot = ACE.Rotations.basis(CA, ll)
   @info("   ... there are $(size(Brot, 2)) rotation-invariance basis functions")

   @info("Find the RPI basis")
   I1 = findall(b -> b[end] == ll, basis.bgrps[1])
   @info("   ... there are $(length(I1)) basis groups with this `ll`:")
   println("grp-idx        kk              ll         length    B-idx")
   for i in I1
      b = basis.bgrps[1][i]
      len = basis.firstb[1][i+1] - basis.firstb[1][i]
      bidx = (basis.firstb[1][i]+1):basis.firstb[1][i+1]
      println(" $(i)   |  $(Int64.(b[2]))  |  $(Int64.(b[3]))  |  $(len) | $(bidx)")
   end
   return nothing
end

"""
create a sub-basis for easier evaluation
"""
struct SubBasis{TB} <: IPBasis
   B::TB
   Ib::Vector{Int}
end

SubBasis(B, Ib::AbstractVector) =
   SubBasis(B, convert(Vector{Int}, collect(Ib)))

(B::SubBasis)(Rs) =
   evaluate(B.B, Rs, zeros(Int16,length(Rs)), 0)[B.Ib]

end


text:
using LinearAlgebra
using YAML

parameters = YAML.load(open("/Users/titou/Documents/PFE/RL-STDP/Julia/Modules/cfg.yml"))

mutable struct Network

    n_exc::Int64
    n_neurons::Int64

    v::Array{Float64}
    u::Array{Float64}
    d::Array{Float64}
    a::Array{Float64}

    connection::Array{Int64,2}
    connection_ee::Array{Int64,2}
    connection_ei::Array{Int64,2}
    connection_ie::Array{Int64,2}
    s_ee::Array{Float64,2}
    s_ei::Array{Float64,2}
    s_ie::Array{Float64,2}
    sd_ee::Array{Float64,2}
    sd_ei::Array{Float64,2}
    trace_e::Array{Float64}
    trace_i::Array{Float64}
    ho::Array{Float64}

    exc::Array{Int64}
    inh::Array{Int64}

    I::Array{Float64}
    da::Float64
    arch::Array{Int64}
    param::Dict

    function connect(n_neurons::Int64, n_exc::Int64, connectivity::Float64,
                     arch::Array{Int64}, inh_per::Float64 ) # parameter
        if arch == [0]
            connection = (rand(n_neurons,n_neurons) .< connectivity)
            n_inh = n_neurons - n_exc
            connection[n_exc+1:n_neurons,n_exc+1:n_neurons] .= zeros(n_inh, n_inh)
        else
            architr = [0; arch]
            connection = zeros(n_neurons,n_neurons)
            for i in 1:length(architr)-2
                si1 = sum(architr[1:i])
                si2 = sum(architr[1:i+1])
                si3 = sum(architr[1:i+2])
                ai2 = architr[i+1]
                ai3 = architr[i+2]
                connection[1+si1:si2,1+si2:si3] .= ones(ai2,ai3)
            end
            archinh = [0;arch[2:end-1]]
            for i in 1:length(archinh)-1
                sii1 = sum(archinh[1:i])
                sii2 = sum(archinh[1:i+1])
                si2 = sum(architr[1:i+1])
                si3 = sum(architr[1:i+2])
                a2 = arch[i+1]
                n_con_inh = floor(Int64,a2*inh_per)
                #mat_ie = ones(a2,a2) .- I(a2)
                mat_ie = zeros(a2,a2)
                for lin in 1:a2
                    ind = rand(1:a2,n_con_inh)
                    ind_f = filter(x -> x != lin, ind)
                    mat_ie[lin,ind_f] .+= 1
                end
                connection[n_exc+sii1+1:n_exc+sii2, si2+1:si3] .= mat_ie
                connection[si2+1:si3, n_exc+sii1+1:n_exc+sii2] .= I(a2)
            end
        end
        exc = collect(1:n_exc)
        inh = collect(n_exc+1:n_neurons)
        return connection, exc, inh
    end

    function Network(arch::Array{Int64}, param::Dict = parameters)
        n_exc = sum(arch)
        n_neurons = n_exc+sum(arch[2:end-1])
        v = param["c"]*ones(n_neurons)
        u = 0.2 * v
        d = [(i<=n_exc) ? param["de"] : param["di"] for i in 1:n_neurons]
        a = [(i<=n_exc) ? param["ae"] : param["ai"] for i in 1:n_neurons]
        n_inh = n_neurons-n_exc
        connection, exc, inh = connect(n_neurons, n_exc, 0.1, arch, param["i_per"])
        connection_ee = connection[exc,exc]
        connection_ei = connection[exc,inh]
        connection_ie = connection[inh,exc]
        # maybe set the Exc-> Inh conn weights to high value
        s_ee = param["wee"] .* ones(n_exc, n_exc)
        s_ei = param["wei"] .* ones(n_exc, n_inh)
        s_ie = param["wie"] .* ones(n_inh, n_exc)
        s_ee .*= connection_ee
        s_ei .*= connection_ei
        s_ie .*= connection_ie
        sd_ee = zeros(n_exc,n_exc)
        sd_ei = zeros(n_exc,n_inh)
        trace_e = zeros(n_exc)
        trace_i = zeros(n_inh)
        ho = zeros(n_neurons)
        #I = 13 .* (rand(n_neurons) .- 0.5)
        I = zeros(n_neurons)
        da = 0.0
        arch = arch
        new(n_exc,n_neurons,v,u,d,a,connection,connection_ee,connection_ei,connection_ie,
            s_ee,s_ei,s_ie,sd_ee,sd_ei,trace_e,trace_i,ho,exc,inh,I,da,arch,param)
    end

    function Network(n_neurons::Int64, connectivity::Float64, n_exc::Int64, param::Dict = parameters)
        v = param["c"]*ones(n_neurons)
        u = 0.2 * v
        d = [(i<=n_exc) ? param["de"] : param["di"] for i in 1:n_neurons]
        a = [(i<=n_exc) ? param["ae"] : param["ai"] for i in 1:n_neurons]
        n_inh = n_neurons-n_exc
        connection, exc, inh = connect(n_neurons, n_exc, connectivity, [0])
        connection_ee = connection[exc,exc]
        connection_ei = connection[exc,inh]
        connection_ie = connection[inh,exc]
        # maybe set the Exc-> Inh conn weights to high value
        s_ee = param["wee"] .* ones(n_exc, n_exc)
        s_ei = param["wei"] .* ones(n_exc, n_inh)
        s_ie = param["wie"] .* ones(n_inh, n_exc)
        s_ee .*= connection_ee
        s_ei .*= connection_ei
        s_ie .*= connection_ie
        sd_ee = zeros(n_exc,n_exc)
        sd_ei = zeros(n_exc,n_inh)
        trace_e = zeros(n_exc)
        trace_i = zeros(n_inh)
        ho = zeros(n_neurons)
        #I = 13 .* (rand(n_neurons) .- 0.5)
        I = zeros(n_neurons)
        da = 0.0
        arch = [0]
        new(n_exc,n_neurons,v,u,d,a,connection,connection_ee,connection_ei,connection_ie,
            s_ee,s_ei,s_ie,sd_ee,sd_ei,trace_e,trace_i,ho,exc,inh,I,da,arch, param)
    end
end

# %% functions

# direct spike of the neuron
function input!(net::Network, input_spikes::Array{Int64})
    net.v[input_spikes] .= 40.0
end

# superthreshold current processed by the neuron
function input!(net::Network, input_spikes::Array{Tuple{Int64,Float64}})
    #net.I = 13 .* (rand(net.n_neurons) .- 0.5)
    net.I = zeros(net.n_neurons)
    max_current = net.param["imax"]
    min_current = net.param["imin"]
    indices = [input_spikes[i][1] for i in eachindex(input_spikes)]
    intensity = [input_spikes[i][2] for i in eachindex(input_spikes)]
    net.I[indices] .+= (max_current-min_current) .* intensity .+ min_current
    net.v[indices] .= @. net.v[indices]+0.5*((0.04*net.v[indices]+5)*net.v[indices]+140-net.u[indices]+net.I[indices])
    net.v[indices] .= @. net.v[indices]+0.5*((0.04*net.v[indices]+5)*net.v[indices]+140-net.u[indices]+net.I[indices])
end

function spike!(net::Network)
    # spikes
    #net.I = 13 .* (rand(net.n_neurons) .- 0.5)
    net.I = zeros(net.n_neurons)
    spiked = findall(x->x>30,net.v .- net.ho)
    spiked_e = filter(x->(x in net.exc), spiked)
    spiked_i = filter(x->(x in net.inh), spiked)
    spiked_ho = filter(x-> x > net.arch[1], spiked_e)
    net.ho[spiked_ho] .+= net.param["theta"]
    net.v[spiked] .= net.param["c"]
    net.u[spiked] .+= net.d[spiked]
    # increment current
    for neuron in spiked
        if neuron in net.exc
            net.I[net.exc] .+= net.s_ee[neuron,:]
            net.I[net.inh] .+= net.s_ei[neuron,:]
        elseif neuron in net.inh
            net.I[net.exc] .+= net.s_ie[neuron-net.n_exc,:]
        end
    end
    #time step (Izhikevic model)
    net.v .= @. net.v+0.5*((0.04*net.v+5)*net.v+140-net.u+net.I)
    net.v .= @. net.v+0.5*((0.04*net.v+5)*net.v+140-net.u+net.I)
    net.u .= @. net.u+net.a*(0.2*net.v-net.u)
    net.trace_e[spiked_e] .= 0.1
    net.trace_i[spiked_i .- net.n_exc] .= 0.1
    # LTP & LTD
    for neuron in spiked_e
        net.sd_ee[:,neuron] .+= net.param["apost"] .* net.trace_e #parameter (a_post)
        net.sd_ee[neuron,:] .-= net.param["apre"] .* net.trace_e # parameter (a_pre)
        if net.arch == [0]
            net.sd_ei[neuron,:] .-= net.param["apre"] .* net.trace_i
        end
    end
    if net.arch == [0]
        for neuron in spiked_i
            net.sd_ei[:,neuron - net.n_exc] .+= net.param["apost"] .* net.trace_e
        end
    end

    # time decay
    net.ho .*= net.param["ttheta"]
    net.trace_e .*= 0.95
    net.trace_i .*= 0.95
    net.da *= net.param["tda"]
    spiked
end

function learn!(net::Network)
    #increase the synaptic weigths
    net.s_ee .+= (0.002+net.da) .* net.sd_ee
    net.s_ee .= clamp.(net.s_ee, 0, net.param["wmax"])
    net.s_ee .*= net.connection_ee
    net.sd_ee .*= 0.99
    if net.arch == [0]
        net.s_ei .+= (0.002+net.da) .* net.sd_ei
        net.s_ei .= clamp.(net.s_ei, 0, net.param["wmax"])
        net.s_ei .*= net.connection_ei
        net.sd_ei .*= 0.99
    end
end

# %% main functions

function step!(net::Network, train::Bool=false)
    spiked = spike!(net)
    if train
        learn!(net)
    end
    spiked
end

function step!(net::Network, input_spikes::Array{Int64}, train::Bool=false)
    input!(net,input_spikes)
    spiked = spike!(net)
    if train
        learn!(net)
    end
    spiked
end

function step!(net::Network, input_spikes::Array{Tuple{Int64,Float64}}, train::Bool=false)
    input!(net,input_spikes)
    spiked = spike!(net)
    if train
        learn!(net)
    end
    spiked
end

# %% visualization function

function weight_viz(net::Network, method::String = "histogram")
    weights = net.s_ee
    if method == "histogram"
        filtered = filter(x-> x!=0.0, weights)
        display(histogram(filtered,
                        legend = false,
                        xlabel = "Synaptic strength",
                        ylabel = "Synapse number"))
    elseif method == "heatmap"
        heatmap(weights,
                title = "connection matrix")
    end
end

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
using HTTP
using JSON
import ARFFFiles
import ScientificTypes: Continuous, Count, Textual, Multiclass, coerce, autotype
using Markdown

const API_URL = "https://www.openml.org/api/v1/json"

# Data API
# The structures are based on these descriptions
# https://github.com/openml/OpenML/tree/master/openml_OS/views/pages/api_new/v1/xsd
# https://www.openml.org/api_docs#!/data/get_data_id

# TODO:
# - Use e.g. DataDeps to cache data locally
# - Put the ARFF parser to a separate package or use ARFFFiles when
#   https://github.com/cjdoris/ARFFFiles.jl/issues/4 is fixed.

"""
Returns information about a dataset. The information includes the name,
information about the creator, URL to download it and more.

- 110 - Please provide data_id.
- 111 - Unknown dataset. Data set description with data_id was not found in the database.
- 112 - No access granted. This dataset is not shared with you.
"""
function load_Dataset_Description(id::Int; api_key::String="")
    url = string(API_URL, "/data/$id")
    try
        r = HTTP.request("GET", url)
        if r.status == 200
            return JSON.parse(String(r.body))
        elseif r.status == 110
            println("Please provide data_id.")
        elseif r.status == 111
            println("Unknown dataset. Data set description with data_id was not found in the database.")
        elseif r.status == 112
            println("No access granted. This dataset is not shared with you.")
        end
    catch e
        println("Error occurred : $e")
        return nothing
    end
    return nothing
end

"""
    MLJOpenML.load(id; parser = :arff)

Load the OpenML dataset with specified `id`, from those listed by
[`list_datasets`](@ref) or on the [OpenML site](https://www.openml.org/search?type=data).
With `parser = :arff` (default) the ARFFFiles.jl parser is used.
With `parser = :auto` the output of the ARFFFiles parser is coerced to
automatically detected scientific types.

Returns a table.

# Examples

```julia
using DataFrames
table = MLJOpenML.load(61);
df = DataFrame(table);
```
"""
function load(id::Int; parser = :arff)
    response = load_Dataset_Description(id)
    arff_file = HTTP.request("GET", response["data_set_description"]["url"])
    data = ARFFFiles.load(IOBuffer(arff_file.body))
    if parser == :auto
        return coerce(data, autotype(data))
    else
        return data
    end
end

"""
Returns a list of all data qualities in the system.

- 412 - Precondition failed. An error code and message are returned
- 370 - No data qualities available. There are no data qualities in the system.
"""
function load_Data_Qualities_List()
    url = string(API_URL, "/data/qualities/list")
    try
        r = HTTP.request("GET", url)
        if r.status == 200
            return JSON.parse(String(r.body))
        elseif r.status == 370
            println("No data qualities available. There are no data qualities in the system.")
        end
    catch e
        println("Error occurred : $e")
        return nothing
    end
    return nothing
end

"""
Returns a list of all data qualities in the system.

- 271 - Unknown dataset. Data set with the given data ID was not found (or is not shared with you).
- 272 - No features found. The dataset did not contain any features, or we could not extract them.
- 273 - Dataset not processed yet. The dataset was not processed yet, features are not yet available. Please wait for a few minutes.
- 274 - Dataset processed with error. The feature extractor has run into an error while processing the dataset. Please check whether it is a valid supported file. If so, please contact the API admins.
"""
function load_Data_Features(id::Int; api_key::String = "")
    if api_key == ""
        url = string(API_URL, "/data/features/$id")
    end
    try
        r = HTTP.request("GET", url)
        if r.status == 200
            return JSON.parse(String(r.body))
        elseif r.status == 271
            println("Unknown dataset. Data set with the given data ID was not found (or is not shared with you).")
        elseif r.status == 272
            println("No features found. The dataset did not contain any features, or we could not extract them.")
        elseif r.status == 273
            println("Dataset not processed yet. The dataset was not processed yet, features are not yet available. Please wait for a few minutes.")
        elseif r.status == 274
            println("Dataset processed with error. The feature extractor has run into an error while processing the dataset. Please check whether it is a valid supported file. If so, please contact the API admins.")
        end
    catch e
        println("Error occurred : $e")
        return nothing
    end
    return nothing
end

"""
Returns the qualities of a dataset.

- 360 - Please provide data set ID
- 361 - Unknown dataset. The data set with the given ID was not found in the database, or is not shared with you.
- 362 - No qualities found. The registered dataset did not contain any calculated qualities.
- 363 - Dataset not processed yet. The dataset was not processed yet, no qualities are available. Please wait for a few minutes.
- 364 - Dataset processed with error. The quality calculator has run into an error while processing the dataset. Please check whether it is a valid supported file. If so, contact the support team.
- 365 - Interval start or end illegal. There was a problem with the interval start or end.
"""
function load_Data_Qualities(id::Int; api_key::String = "")
    if api_key == ""
        url = string(API_URL, "/data/qualities/$id")
    end
    try
        r = HTTP.request("GET", url)
        if r.status == 200
            return JSON.parse(String(r.body))
        elseif r.status == 360
            println("Please provide data set ID")
        elseif r.status == 361
            println("Unknown dataset. The data set with the given ID was not found in the database, or is not shared with you.")
        elseif r.status == 362
            println("No qualities found. The registered dataset did not contain any calculated qualities.")
        elseif r.status == 363
            println("Dataset not processed yet. The dataset was not processed yet, no qualities are available. Please wait for a few minutes.")
        elseif r.status == 364
            println("Dataset processed with error. The quality calculator has run into an error while processing the dataset. Please check whether it is a valid supported file. If so, contact the support team.")
        elseif r.status == 365
            println("Interval start or end illegal. There was a problem with the interval start or end.")
        end
    catch e
        println("Error occurred : $e")
        return nothing
    end
    return nothing
end

"""
    load_List_And_Filter(filters; api_key = "")

See [OpenML API](https://www.openml.org/api_docs#!/data/get_data_list_filters).
"""
function load_List_And_Filter(filters::String; api_key::String = "")
    if api_key == ""
        url = string(API_URL, "/data/list/$filters")
    end
    try
        r = HTTP.request("GET", url)
        if r.status == 200
            return JSON.parse(String(r.body))
        elseif r.status == 370
            println("Illegal filter specified.")
        elseif r.status == 371
            println("Filter values/ranges not properly specified.")
        elseif r.status == 372
            println("No results. There where no matches for the given constraints.")
        elseif r.status == 373
            println("Can not specify an offset without a limit.")
        end
    catch e
        println("Error occurred : $e")
        return nothing
    end
    return nothing
end

qualitynames(x) = haskey(x, "name") ? [x["name"]] : []

"""
    list_datasets(; tag = nothing, filters = "" api_key = "", output_format = NamedTuple)

Lists all active OpenML datasets, if `tag = nothing` (default).
To list only datasets with a given tag, choose one of the tags in [`list_tags()`](@ref).
An alternative `output_format` can be chosen, e.g. `DataFrame`, if the
`DataFrames` package is loaded.

A filter is a string of `<data quality>/<range>` or `<data quality>/<value>`
pairs, concatenated using `/`, such as

```julia
    filter = "number_features/10/number_instances/500..10000"
```

The allowed data qualities include `tag`, `status`, `limit`, `offset`,
`data_id`, `data_name`, `data_version`, `uploader`,
`number_instances`, `number_features`, `number_classes`,
`number_missing_values`.

For more on the format and effect of `filters` refer to the [openml
API](https://www.openml.org/api_docs#!/data/get_data_list_filters).

# Examples
```
julia> using DataFrames

julia> ds = MLJOpenML.list_datasets(
               tag = "OpenML100",
               filter = "number_instances/100..1000/number_features/1..10",
               output_format = DataFrame
)

julia> sort!(ds, :NumberOfFeatures)
```
"""
function list_datasets(; tag = nothing, filter = "", filters=filter,
                         api_key = "", output_format = NamedTuple)
    if tag !== nothing
        if is_valid_tag(tag)
            filters *= "/tag/$tag"
        else
            @warn "$tag is not a valid tag. See `list_tags()` for a list of tags."
            return
        end
    end
    data = MLJOpenML.load_List_And_Filter(filters; api_key = api_key)
    datasets = data["data"]["dataset"]
    qualities = Symbol.(union(vcat([vcat(qualitynames.(entry["quality"])...) for entry in datasets]...)))
    result = merge((id = Int[], name = String[], status = String[]),
                   NamedTuple{tuple(qualities...)}(ntuple(i -> Union{Missing, Int}[], length(qualities))))
    for entry in datasets
        push!(result.id, entry["did"])
        push!(result.name, entry["name"])
        push!(result.status, entry["status"])
        for quality in entry["quality"]
            push!(getproperty(result, Symbol(quality["name"])),
                  Meta.parse(quality["value"]))
        end
        for quality in qualities
            if length(getproperty(result, quality)) < length(result.id)
                push!(getproperty(result, quality), missing)
            end
        end
    end
    output_format(result)
end

is_valid_tag(tag::String) = tag ∈ list_tags()
is_valid_tag(tag) = false

"""
    list_tags()

List all available tags.
"""
function list_tags()
    url = string(API_URL, "/data/tag/list")
    try
        r = HTTP.request("GET", url)
        return JSON.parse(String(r.body))["data_tag_list"]["tag"]
    catch
        return nothing
    end
end

"""
    describe_dataset(id)

Load and show the OpenML description of the data set `id`.
Use [`list_datasets`](@ref) to browse available data sets.

# Examples
```
julia> MLJOpenML.describe_dataset(6)
  Author: David J. Slate Source: UCI
  (https://archive.ics.uci.edu/ml/datasets/Letter+Recognition) - 01-01-1991 Please cite: P.
  W. Frey and D. J. Slate. "Letter Recognition Using Holland-style Adaptive Classifiers".
  Machine Learning 6(2), 1991

    1. TITLE:

  Letter Image Recognition Data

  The objective is to identify each of a large number of black-and-white
  rectangular pixel displays as one of the 26 capital letters in the English
  alphabet.  The character images were based on 20 different fonts and each
  letter within these 20 fonts was randomly distorted to produce a file of
  20,000 unique stimuli.  Each stimulus was converted into 16 primitive
  numerical attributes (statistical moments and edge counts) which were then
  scaled to fit into a range of integer values from 0 through 15.  We
  typically train on the first 16000 items and then use the resulting model
  to predict the letter category for the remaining 4000.  See the article
  cited above for more details.
```
"""
describe_dataset(id) =  Markdown.parse(load_Dataset_Description(id)["data_set_description"]["description"])

# Flow API

# Task API

# Run API


text:
abstract type AbstractAsymptoticPosterior{P,T,Pm1,R} <: AbstractDifferentiableObject{P,T} end

struct AsymptoticPosteriorFD{P,T,Pm1,R,D<:AbstractDifferentiableObject{Pm1,T},M,S,MAP <: MAPFD{P,T,R},
                            R2,SM <: MutableFixedSizeMatrix{R2,Pm1,T}} <: AbstractAsymptoticPosterior{P,T,Pm1,R}
    od::D
    method::M
    state::S
    map::MAP
    nlmax::Base.RefValue{T}
    rstar::Base.RefValue{T}
    α::Base.RefValue{T}
    Lsmall::SM
end

## abstract interface
## required methods

@inline set_profile_val_ind!(ap::AsymptoticPosteriorFD, v, i) = set_profile_val_ind!(ap.od.config.f, v, i)
@inline setswap!(ap::AsymptoticPosteriorFD, i) = ap.map.od.config.f.i[] = i
@inline rstar(ap::AsymptoticPosteriorFD) = ap.rstar[]
@inline rstar!(ap::AsymptoticPosteriorFD{P,T}, v::T) where {P,T} = (ap.rstar[] = v)
@inline α(ap::AsymptoticPosteriorFD) = ap.α[]
@inline function α!(ap::AsymptoticPosteriorFD{P,T}, v::T) where {P,T}
    ap.α[] = v
    rstar!(ap, Φ⁻¹(v))
    v
end
@inline mode(ap::AsymptoticPosteriorFD) = ap.map.θhat
@inline mode(ap::AsymptoticPosteriorFD, i::Integer) = @inbounds ap.map.θhat[i]
@inline std_estimates(ap::AsymptoticPosteriorFD) = ap.map.std_estimates
@inline std_estimates(ap::AsymptoticPosteriorFD, i) = @inbounds ap.map.std_estimates[i]

@inline gradient(ap::AsymptoticPosteriorFD) = ap.map.od.config.result.grad
@inline hessian(ap::AsymptoticPosteriorFD) = DiffResults.hessian(ap.map.od.config.result)

@inline profile_ind(ap::AsymptoticPosteriorFD) = ap.od.config.f.i[]
@inline profile_ind!(ap::AsymptoticPosteriorFD, i::Integer) = (ap.od.config.f.i[] = i)
@inline profile_val(ap::AsymptoticPosteriorFD) = ap.od.config.f.v[]
@inline nl_profile_max(ap::AsymptoticPosteriorFD) = ap.nlmax[]
@inline nl_profile_max!(ap::AsymptoticPosteriorFD{P,T}, v::T) where {P,T} = (ap.nlmax[] = v)
@inline nl_max(ap::AsymptoticPosteriorFD) = ap.map.nlmax[]
@inline nl_max!(ap::AsymptoticPosteriorFD{P,T}, v::T) where {P,T} = (ap.map.nlmax[] = v)

@inline fit!(ap::AsymptoticPosteriorFD{P,T}, x::AbstractFixedSizeVector{P,T}) where {P,T} = fit!(ap.map, x)

"""
initial_val_buffer(ap::AbstractAsymptoticPosterior{P,T}) where {P,T}
returns a mutable buffer of length {P}.
The buffer shuold be <: AbstractArray{T}.
Each call must refer to the same piece of memory.
"""
@inline initial_val_buffer(ap::AsymptoticPosteriorFD) = ap.map.buffer

@inline nuisance_parameters(ap::AsymptoticPosteriorFD) = DifferentiableObjects.ref_x_old(ap.state)

@inline base_adjustment(ap::AsymptoticPosteriorFD) = ap.map.base_adjust[]
@inline base_adjustment!(ap::AsymptoticPosteriorFD, v) = (ap.map.base_adjust[] = v)
@inline Lfull(ap::AsymptoticPosteriorFD) = ap.map.Lfull

function refit!(ap::AbstractAsymptoticPosterior{P}) where {P}
    i = profile_ind(ap)
    setswap!(ap, P)
    fit!(ap.map, ap.od.config.f.x)
    setswap!(ap, i)
end

@inline function delta_log_likelihood(ap::AbstractAsymptoticPosterior)#{P}) where {P}
#    return nl_profile_max(ap) - nl_max(ap)
    δ = nl_profile_max(ap) - nl_max(ap)
    if δ < 0
        refit!(ap)
        return delta_log_likelihood(ap)
#        @warn "δ log likelihood = $δ < 0, returning 0."
#        return zero(δ)
    end
    δ
end

@inline profile_lower_triangle(ap::AsymptoticPosteriorFD) = ap.Lsmall

@inline function optimize_profile!(ap::AsymptoticPosteriorFD{P,T}) where {P,T}
    scaled_nlmax, scale = optimize_scale!(ap.state, ap.od, nuisance_parameters(ap), ap.method, one(T))
    nl_profile_max!(ap, scaled_nlmax / scale)
end

@inline function DifferentiableObjects.hessian!(ap::AsymptoticPosteriorFD)
    DifferentiableObjects.hessian!(ap.map.od.config, ap.map.buffer)
end

##
## Abstract interface, optional methods
## These are optional, because they all have the default fallbacks, shown below.
##

"""
Naively assumes indepedence, and just uses the global maximum for the profile maximum.
This seems like it ought to perform poorly, but in tests it appears more robust.
"""
function naive_expected_nuisance!(ap::AbstractAsymptoticPosterior{P}, i::Int = profile_ind(ap)) where {P}
    nuisance = nuisance_parameters(ap)
    θhat = mode(ap)
    @inbounds for j in 1:i-1
        nuisance[j] = θhat[j]# + ap.map.cov[j,i] * sf
    end
    @inbounds for j in i:P-1
        nuisance[j] = θhat[j+1]# + ap.map.cov[j,i] * sf
    end
end

function profilepdf(ap::AbstractAsymptoticPosterior{P,T}, ::Val{reset_search_start}=Val{true}()) where {reset_search_start,P,T}
    i = profile_ind(ap)
    reset_search_start && naive_expected_nuisance!(ap, i)
    optimize_profile!(ap)
end
function profilepdf(ap::AbstractAsymptoticPosterior{P,T}, x, i, ::Val{reset_search_start}=Val{true}()) where {reset_search_start,P,T}
    set_profile_val_ind!(ap, x, i)
    # Reseset the other parameters to the corresponding mode.
    reset_search_start && naive_expected_nuisance!(ap, i)
    optimize_profile!(ap)
end

function (ap::AbstractAsymptoticPosterior)(theta, i = profile_ind(ap))
    debug() && println("")
    debug() && @show theta, rstar(ap), i
    debug() && @assert isfinite(theta)
    debug() && @show rstar_p(ap, theta, i) + rstar(ap)
    rstar_p(ap, theta, i) + rstar(ap)
end

# Positive if x is smaller than mode, negative if bigger
function rp(ap::AsymptoticPosteriorFD, x, i = profile_ind(ap))
    copysign( sqrt(2(delta_log_likelihood(ap))), mode(ap, i) - x)
end

function rstar_p(ap::AbstractAsymptoticPosterior, theta::Number, i::Int=profile_ind(ap))

    profilepdf(ap, theta, i)
    set_buffer_to_profile!(ap, i)

    setswap!(ap, i)
    hessian!(ap)

    debug() && @show ap.map.buffer
    debug() && @show hessian(ap)

    r = rp(ap, theta, i)
    q = qb(ap, theta, i)
#    @show r, q
    r + log(q/r)/r
end

sym(a,i) = Symbol(a, :_, i)
"""
Generates quote for vectorized multiplication of
    (L * g)' * (L * h)

"""
function profile_correction_quote(P, R, T)
    pad = R - P + 1
    # @show P,R
    VL = VectorizationBase.pick_vector_width(R, T)
    # VL = min(R, VectorizationBase.REGISTER_SIZE ÷ sizeof(T))
    VLT = VL * sizeof(T)
    V = Vec{VL,T}
    # q = quote @fastmath @inbounds begin end end
    # qa = q.args[2].args[3].args[3].args
    q = quote $(Expr(:meta, :inline)) end
    qa = q.args
    iter = R ÷ VL # number of iterations down the columns of the matrix.
    push!(qa, :(ptr_Li = pointer(Li)))
    push!(qa, :(@inbounds $(sym(:vH,0)) = SIMDPirates.vbroadcast($V,hess[1,$P]) ))
    push!(qa, :(@inbounds $(sym(:vG,0)) = SIMDPirates.vbroadcast($V,grad[1]) ))
    # r = 0
    push!(qa, :($(sym(:vL,0)) = vload($V, ptr_Li) ))
    push!(qa, :($(sym(:vHL,0)) = SIMDPirates.vmul( $(sym(:vH,0)), $(sym(:vL,0)) )))
    push!(qa, :($(sym(:vGL,0)) = SIMDPirates.vmul($(sym(:vG,0)), $(sym(:vL,0)) )))
    for r ∈ 1:iter-1 #ps = 0
        push!(qa, :($(sym(:vL,r)) = vload($V, ptr_Li + $(r*VLT)) ))
        push!(qa, :($(sym(:vHL,r)) = SIMDPirates.vmul($(sym(:vH,0)),$(sym(:vL,r)) )))
        push!(qa, :($(sym(:vGL,r)) = SIMDPirates.vmul($(sym(:vG,0)),$(sym(:vL,r)) )))
    end
    for ps ∈ 1:VL-pad-1
        push!(qa, :(@inbounds $(sym(:vH,ps)) = SIMDPirates.vbroadcast($V,hess[$(ps+1),$P])))
        push!(qa, :(@inbounds $(sym(:vG,ps)) = SIMDPirates.vbroadcast($V,grad[$(ps+1)])))
        for r ∈ 0:iter-1
            push!(qa, :($(sym(:vL,r)) = vload($V, ptr_Li + $(r*VLT + sizeof(T)*R*ps))))
            push!(qa, :($(sym(:vHL,r)) = SIMDPirates.vmuladd($(sym(:vH,ps)),$(sym(:vL,r)),$(sym(:vHL,r))) ))
            push!(qa, :($(sym(:vGL,r)) = SIMDPirates.vmuladd($(sym(:vG,ps)),$(sym(:vL,r)),$(sym(:vGL,r))) ))
        end
    end
    push!(qa, :(Vout = SIMDPirates.vmul($(sym(:vHL,0)), $(sym(:vGL,0)))))
    itermin = 1
    # for pb ∈ 0:VL:P-2
    #     for ps ∈ max(1,pb):min( pb+VL-1, P-2 )
    for pb ∈ VL-pad:VL:P-2
        for ps ∈ pb:min( pb+VL-1, P-2 )
            push!(qa, :(@inbounds $(sym(:vH,ps)) = SIMDPirates.vbroadcast($V,hess[$(ps+1),$P])))
            push!(qa, :(@inbounds $(sym(:vG,ps)) = SIMDPirates.vbroadcast($V,grad[$(ps+1)])))
            for r ∈ itermin:iter-1
                push!(qa, :($(sym(:vL,r)) = vload($V, ptr_Li + $(r*VLT + sizeof(T)*R*ps))))
                push!(qa, :($(sym(:vHL,r)) = SIMDPirates.vmuladd($(sym(:vH,ps)),$(sym(:vL,r)),$(sym(:vHL,r))) ))
                push!(qa, :($(sym(:vGL,r)) = SIMDPirates.vmuladd($(sym(:vG,ps)),$(sym(:vL,r)),$(sym(:vGL,r))) ))
            end
        end
        push!(qa, :(Vout = SIMDPirates.vmuladd($(sym(:vHL,itermin)), $(sym(:vGL,itermin)),Vout)))
        itermin += 1
    end
    # push!(qa, :(println("Finished profile correction.")))
    push!(qa, :(SIMDPirates.vsum(Vout)))
    q
end
"""
Calculates the quadratic form of grad' * Li' * Li * hess[1:end-1,end]
"""
@generated function profile_correction(Li::PaddedMatrices.AbstractMutableFixedSizeMatrix{R,Pm1,T,R},
                        grad::PaddedMatrices.AbstractFixedSizeVector, hess::PaddedMatrices.AbstractFixedSizeMatrix) where {Pm1,T,R}
    profile_correction_quote(Pm1+1, R, T)
end

function fdf_adjrstar_p(
    ap::AbstractAsymptoticPosterior{P,T}, theta, p_i::Int=profile_ind(ap),
    ::Val{reset_search_start} = Val{true}(), repeats = 0
) where {P,T,reset_search_start}

    profilepdf(ap, theta, p_i, Val{reset_search_start}())
    set_buffer_to_profile!(ap, p_i)
    setswap!(ap, p_i)
    hess = hessian!(ap)

    delta_log_like = delta_log_likelihood(ap)
    @inbounds r = copysign( sqrt((T(2))*delta_log_like), mode(ap, p_i) - theta)
    # r = rp(ap, theta, p_i)

    grad = gradient(ap)

    Li = profile_lower_triangle(ap)
    rootdet = PaddedMatrices.invcholdetLLc!(Li, hess)

    prof_factor = profile_correction(Li, grad, hess)
    hess_adjust = rootdet * base_adjustment(ap)
    @inbounds q = (prof_factor - grad[P]) * hess_adjust
    #    @show r, q, delta_log_like, rootdet, base_adjustment(ap), hess_adjust
    if sign(q)  != sign(r)
        repeats > 10 && return T(-Inf)
        refit!(ap)
        return fdf_adjrstar_p(ap, theta, p_i, Val{reset_search_start}(), repeats + 1)
    end
    r⭐ = r + log(q/r)/r
#    r⭐ + rstar(ap), exp(T(0.5)*abs2(r⭐)-delta_log_like) / hess_adjust#, exp(-delta_log_like) / hess_adjust
    1 - Φ(r⭐) - α(ap), exp(-delta_log_like) / (Base.FastMath.sqrt_fast(2π) * rootdet * base_adjustment(ap))
#    Φ(r⭐ + rstar(ap)), exp(delta_log_like) / (Base.FastMath.sqrt_fast(2π) * rootdet * base_adjustment(ap))
end

function cdf(
    ap::AbstractAsymptoticPosterior{P,T}, theta, p_i::Int=profile_ind(ap),
    ::Val{reset_search_start} = Val{true}(), tries = 0
) where {P,T,reset_search_start}

    profilepdf(ap, theta, p_i, Val{reset_search_start}())
    set_buffer_to_profile!(ap, p_i)
    setswap!(ap, p_i)
    hess = hessian!(ap)
    delta_log_like = delta_log_likelihood(ap)
    @inbounds r = copysign( sqrt((T(2))*delta_log_like), mode(ap, p_i) - theta)
    grad = gradient(ap)
    Li = profile_lower_triangle(ap)
    rootdet = PaddedMatrices.invcholdetLLc!(Li, hess)
    prof_factor = profile_correction(Li, grad, hess)
    hess_adjust = rootdet * base_adjustment(ap)
    @inbounds q = (prof_factor - grad[P]) * hess_adjust
    if sign(q) != sign(r)
        refit!(ap)
        if tries > 10
            @warn "Tried refitting 10 times."
            return -Inf
        end
        return cdf(ap, theta, p_i, Val{reset_search_start}(), tries + 1)
    end
    r⭐ = r + log(q/r)/r

    1 - Φ(r⭐)
end
@generated function subhessian(H::AbstractMutableFixedSizeMatrix{P,P,T,R}) where {P,T,R}
    Pm1 = P - 1
    quote
        $(Expr(:meta,:inline))
        PtrMatrix{$Pm1,$Pm1,$T,$R,$(Pm1*R)}(pointer(H))
    end
end

function pdf(ap::AbstractAsymptoticPosterior{P,T,Pm1}, theta, i = profile_ind(ap),::Val{reset_search_start}=Val{true}()) where {P,T,Pm1,reset_search_start}

    profilepdf(ap, theta, i, Val{reset_search_start}())
    set_buffer_to_profile!(ap, i)
    setswap!(ap, i)
    hessian!(ap.map.od.config, ap.map.buffer)

    # is there a more efficient way of calculating the determinant?
    # skip storing? Implement that...
    # rootdet = PaddedMatrices.choldet!(hessian(ap), Val{Pm1}())
    # isfinite(rootdet) ? exp(delta_log_likelihood(ap)) / (sqrt(2π) * rootdet * base_adjustment(ap)) : T(Inf)
    rootdet, success = PaddedMatrices.safecholdet!(subhessian(hessian(ap)))
    success ? exp(-delta_log_likelihood(ap)) / (sqrt(2π) * rootdet * base_adjustment(ap)) : T(-Inf)
end

function qb(ap::AbstractAsymptoticPosterior{P}, theta, pi = profile_ind(ap)) where {P}
    grad = gradient(ap)
    hess = hessian(ap) # P x P
    Li = profile_lower_triangle(ap)
    rootdet = PaddedMatrices.invcholdetLLc!(Li, hess)
    prof_factor = profile_correction(Li, grad, hess)

    @inbounds (prof_factor - grad[P]) * rootdet * base_adjustment(ap)
end

function set_buffer_to_profile!(ap::AbstractAsymptoticPosterior{P}, i = profile_ind(ap)) where {P}
    nuisance = nuisance_parameters(ap)
    buffer = initial_val_buffer(ap)
    @inbounds for j in 1:i-1
        buffer[j] = nuisance[j]
    end
    if i != P
        @inbounds buffer[i] = nuisance[P-1]
        @inbounds for j in i+1:P-1
            buffer[j] = nuisance[j-1]
        end
    end
    buffer[P] = profile_val(ap)
    nothing
end

Φ(x::T) where {T} = muladd(T(0.5), erf( x / Base.FastMath.sqrt_fast(2) ), T(0.5))
Φ⁻¹(x::T) where {T} = Base.FastMath.sqrt_fast(T(2))*erfinv(muladd(T(2),x,T(-1)))
function Statistics.quantile(ap::AbstractAsymptoticPosterior, alpha, i)
    profile_ind!(ap, i)
    α!(ap, alpha)
#    rstar!(ap, Φ⁻¹(alpha))
    # quadratic_search(ap, i)
    linear_search(ap, profile_ind(ap))
end
function Statistics.quantile(ap::AbstractAsymptoticPosterior, alpha)
    α!(ap, alpha)
#    rstar!(ap, Φ⁻¹(alpha))
    linear_search(ap, profile_ind(ap))
    # quadratic_search(ap, profile_ind(ap))
end

export lquantile, qquantile
function lquantile(ap::AbstractAsymptoticPosterior, alpha, i)
    profile_ind!(ap, i)
    α!(ap, alpha)
    linear_search(ap, i)
end
function lquantile(ap::AbstractAsymptoticPosterior, alpha)
    α!(ap, alpha)
#    rstar!(ap, Φ⁻¹(alpha))
    linear_search(ap, profile_ind(ap))
end

function qquantile(ap::AbstractAsymptoticPosterior, alpha, i)
    profile_ind!(ap, i)
    α!(ap, alpha)
#    rstar!(ap, Φ⁻¹(alpha))
    quadratic_search(ap, i)
end
function qquantile(ap::AbstractAsymptoticPosterior, alpha)
    α!(ap, alpha)
#    rstar!(ap, Φ⁻¹(alpha))
    quadratic_search(ap, profile_ind(ap))
end

# (ap::AsymptoticPosteriorFD)(x) = ap.map.od(x)


text:
"""
Fast, lightweight quaternion module.
"""
module Quaternions

using StaticArrays: SVector, MVector
using LinearAlgebra

export Quaternion, FreeQuat, UnitQuat

######################## Quat #############################

abstract type NormConstraint end
struct FreeNorm <: NormConstraint end
struct UnitNorm <: NormConstraint end

struct Quaternion{N <: NormConstraint} <: AbstractVector{Float64}
    _sv::SVector{4, Float64}
    function Quaternion{N}(v::AbstractVector, normalization::Bool) where {N}
        sv = SVector{4,Float64}(v) #faster normalization method
        return normalization ? new{N}(normalize(sv)) : new{N}(sv)
    end
end

const FreeQuat = Quaternion{FreeNorm}
const UnitQuat = Quaternion{UnitNorm}

Quaternion{FreeNorm}(v::AbstractVector; normalization::Bool = false) = FreeQuat(v, normalization)
Quaternion{UnitNorm}(v::AbstractVector; normalization::Bool = true) = UnitQuat(v, normalization)

Quaternion{N}(s::Real) where {N} = convert(Quaternion{N}, s)

function Quaternion{FreeNorm}(; real = 0.0, imag = zeros(SVector{3}))
    FreeQuat(vcat(real, SVector{3,Float64}(imag)))
end

function Quaternion{UnitNorm}(; real = 0.0, imag = zeros(SVector{3}), normalization::Bool = true)
    sv = vcat(real, SVector{3,Float64}(imag))
    any(sv .!= 0) || throw(ArgumentError("UnitQuat requires a non-zero real or imaginary part"))
    UnitQuat(vcat(real, SVector{3,Float64}(imag)); normalization)
end

#### AbstractVector interface
Base.size(::Quaternion) = (4,)
Base.length(::Quaternion) = 4
Base.firstindex(::Quaternion) = 1
Base.lastindex(::Quaternion) = 4
Base.eltype(::Quaternion) = Float64

#### Indexing and retrieving real and imaginary parts
Base.getindex(q::Quaternion, i) = getfield(q, :_sv)[i]
Base.getindex(q::Quaternion, s::Symbol) = getindex(q, Val(s))
Base.getindex(q::Quaternion, ::Val{:real}) = getfield(q, :_sv)[1]
Base.getindex(q::Quaternion, ::Val{:imag}) = SVector{3, Float64}(@view getfield(q, :_sv)[2:4])
Base.getproperty(q::Quaternion, s::Symbol) = getindex(q, Val(s))

Base.copy(q::Quaternion{N}) where {N} = Quaternion{N}(copy(getfield(q, :_sv)))
Base.show(io::IO, q::Quaternion) = print(io, "$(typeof(q))($(q[:]))")
Base.show(io::IO, ::MIME"text/plain", q::Quaternion) = print(io, "$(typeof(q))($(q[:]))")

#### Normalization
LinearAlgebra.norm(q::Quaternion) = norm(getfield(q, :_sv)) #uses StaticArrays implementation
LinearAlgebra.normalize(f::FreeQuat) = FreeQuat(normalize(getfield(f, :_sv)))
LinearAlgebra.normalize(u::UnitQuat) = UnitQuat(getfield(u, :_sv))
@inline norm_sqr(q::Quaternion) = (data = getfield(q,:_sv); sum(data.*data))

#### Conjugate & Inverse
@inline Base.adjoint(q::Quaternion) = conj(q)

@inline Base.conj(q::FreeQuat) = FreeQuat(vcat(q.real, -q.imag))
@inline Base.inv(q::FreeQuat) = FreeQuat(getfield(q', :_sv) / norm_sqr(q))

@inline Base.conj(u::UnitQuat) = UnitQuat(vcat(u.real, -u.imag), false)
@inline Base.inv(u::UnitQuat) = u'

#### Conversion & Promotion
Base.convert(::Type{Quaternion{N}}, s::Real) where {N} = Quaternion{N}(SVector{4,Float64}(s, 0, 0, 0))
Base.convert(::Type{Quaternion{N}}, v::AbstractVector) where {N} = Quaternion{N}(v)

Base.promote_rule(::Type{<:Quaternion}, ::Type{S}) where {S<:Real} = FreeQuat
Base.promote_rule(::Type{Quaternion{N}}, ::Type{<:AbstractVector{<:Real}}) where {N} = Quaternion{N}
Base.promote_rule(::Type{UnitQuat}, ::Type{FreeQuat}) = FreeQuat

#### Operators
Base.:+(q::FreeQuat) = q
Base.:-(q::FreeQuat) = FreeQuat(-getfield(q, :_sv))
Base.:+(u::UnitQuat) = u
Base.:-(u::UnitQuat) = UnitQuat(-getfield(u, :_sv), false)

Base.:(==)(q1::Quaternion, q2::Quaternion) = getfield(q1,:_sv) == getfield(q2,:_sv)
Base.:(≈)(q1::Quaternion, q2::Quaternion) = getfield(q1,:_sv) ≈ getfield(q2,:_sv)

Base.:+(q1::FreeQuat, q2::FreeQuat) = FreeQuat(getfield(q1,:_sv) + getfield(q2,:_sv))
Base.:+(q1::Quaternion, q2::Quaternion) = +(promote(q1,q2)...)
Base.:+(q::Quaternion, a::Real) = +(promote(q, a)...)
Base.:+(a::Real, q::Quaternion) = +(promote(a, q)...)

Base.:-(q1::FreeQuat, q2::FreeQuat) = FreeQuat(getfield(q1,:_sv) - getfield(q2,:_sv))
Base.:-(q1::Quaternion, q2::Quaternion) = -(promote(q1,q2)...)
Base.:-(q::Quaternion, a::Real) = -(promote(q, a)...)
Base.:-(a::Real, q::Quaternion) = -(promote(a, q)...)

Base.:*(q1::Quaternion{N1}, q2::Quaternion{N2}) where {N1,N2} = *(promote(q1,q2)...)

@inline function Base.:*(q1::Quaternion{N}, q2::Quaternion{N}) where {N}

    p_re = q1.real * q2.real - q1.imag ⋅ q2.imag
    p_im = q1.real * q2.imag + q2.real * q1.imag + q1.imag × q2.imag

    Quaternion{N}(vcat(p_re, p_im), false) #bypass normalization for UnitQuat!
end

Base.:*(q::Quaternion, a::Real) = a * q
Base.:*(a::Real, q::Quaternion) = FreeQuat(a * getfield(q, :_sv))

Base.:/(q1::Quaternion{N}, q2::Quaternion{N}) where {N} = Quaternion{N}(q1 * inv(q2))
Base.:/(q1::Quaternion{N1}, q2::Quaternion{N2}) where {N1,N2} = /(promote(q1,q2)...)
Base.:/(q::Quaternion, a::Real) = FreeQuat(getfield(q, :_sv)/ a)
Base.:/(a::Real, q::Quaternion) = /(promote(a, q)...)

Base.:\(q1::Quaternion{N}, q2::Quaternion{N}) where {N} = Quaternion{N}(inv(q1) * q2) #!= /(q2, q1) == q2 * inv(q1)
Base.:\(q1::Quaternion{N1}, q2::Quaternion{N2}) where {N1,N2} = \(promote(q1,q2)...)
Base.:\(q::Quaternion, a::Real) = \(promote(q, a)...)
Base.:\(a::Real, q::Quaternion) = q / a

end


text:
abstract type BooleanGenerator <: DiscreteFunctionGenerator end
const BoolGen = BooleanGenerator

"""
Generates all boolean function of arity = n
"""
struct BooleanFunctionGenerator{T, F} <: BoolGen
    itr :: UnitRange{T}
    farr :: Vector{F}
end

struct BooleanMapGenerator{T, F} <: BoolGen
    itr :: T
    farr :: Matrix{F}
end

const BFGen = BooleanFunctionGenerator
const BMGen = BooleanMapGenerator
# from arity
function BFGen(arity)
    # we want to iterate over all functions
    itr = 0 : (2^(2^arity) - 1)
    farr = zeros(Int, 2^arity)
    return BFGen(itr, farr)
end

function BMGen(arity, nfuns)
    # we want to iterate over all functions
    fun_itr = 0 : (2^(2^arity) - 1)
    itr = product(fill(fun_itr, nfuns)...)
    farr = zeros(Int, nfuns, 2^arity)
    return BMGen(itr, farr)
end

BFGen(bc :: BCube) = BFGen(ndims(bc))

nfuns(gen :: BFGen) = 1
nfuns(gen :: BMGen) = size(gen.farr)[1]
arity(gen :: BFGen) = Int(log2(length(gen.farr)))
arity(gen :: BFGen) = Int(log2(size(gen.farr)[2]))

# iteration utilities #
function onestep(ctr, gen :: BFGen)
    value = ctr[1]
    state = ctr[2]
    fill_bits!(gen.farr, value)
    f = BoolFun(gen.farr)
    return (f, state)
end

function onestep(ctr, gen :: BMGen)
    value = ctr[1]
    state = ctr[2]
    for i in 1 : nfuns(gen)
        fill_bits!(view(gen.farr, i, :), value[i])
    end
    f = BoolMap(gen.farr)
    return (f, state)
end

onestep(:: Nothing, gen :: BFGen) = nothing
onestep(:: Nothing, gen :: BMGen) = nothing
Base.iterate(gen :: BoolGen) = onestep(iterate(gen.itr), gen)
Base.iterate(gen :: BoolGen, state) = onestep(iterate(gen.itr, state), gen)

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions