text:
# Macros

"""
    autohasheq(typ)

Generate hash and equality methods for a type that are sensitive to type parameters and field
values. This macro is a modified version of AutoHashEquals.jl, where type parameters are not taken
into account.
"""
macro autohasheq(typ)
    @assert typ.head == :struct

    typename = nothing
    typeargs = []
    fieldnames = Any[]

    # Deconstruct the type declaration to get name and type parameters
    declaration = typ.args[2]
    if declaration.head == :<:
        declaration = declaration.args[1]
    end
    if isa(declaration, Symbol)
        typename = declaration
    elseif declaration.head == :curly
        typename = declaration.args[1]
        typeargs = declaration.args[2:end]
    end

    # Deconstruct field names
    for line in typ.args[3].args
        isa(line, LineNumberNode) && continue
        isa(line, Symbol) && (push!(fieldnames, line); continue)
        if line.head == :(::)
            push!(fieldnames, line.args[1])
        end
    end

    # Hash function uses type name, type parameters and field values
    hashcode = :(hash($(QuoteNode(typename)), h))
    for typearg in reverse(typeargs)
        hashcode = :(hash($typearg, $hashcode))
    end
    for fieldname in reverse(fieldnames)
        hashcode = :(hash(self.$fieldname, $hashcode))
    end

    partype = isempty(typeargs) ? typename : :($typename{$(typeargs...)})
    hashcode = :(Base.hash(self::$partype, h::UInt) where {$(typeargs...)} = $hashcode)

    # Equality function uses type parameters and field values
    eqcode = :(true)
    for fieldname in reverse(fieldnames)
        eqcode = :(isequal(l.$fieldname, r.$fieldname) && $eqcode)
    end
    eqcode = :(Base.:(==)(l::$partype, r::$partype) where {$(typeargs...)} = $eqcode)

    ineqcode = isempty(typeargs) ? :() : :(Base.:(==)(l::$typename, r::$typename) = false)

    quote
        $(esc(typ))
        $(esc(hashcode))
        $(esc(eqcode))
        $(esc(ineqcode))
    end
end

function _destructure(prototype)
    funcname = prototype.args[1]
    funcparams = prototype.args[2:end]

    destructure, newparams = Any[], Any[]

    for arg in funcparams
        if isa(arg, Symbol) || arg.head == :(::)
            push!(newparams, esc(arg))
            continue
        end

        @assert arg.head == :vect

        parameter = gensym()
        push!(newparams, esc(parameter))

        for (i, subarg) in enumerate(arg.args)
            if isa(subarg, Symbol)
                push!(destructure, esc(:($subarg = $parameter[$i])))
            else
                @assert subarg.head == :...
                push!(destructure, esc(:($(subarg.args[1]) = $parameter[$i:end])))
            end
        end
    end

    return :($(esc(funcname))($(newparams...))), destructure
end

macro destructure(funcdef)
    @assert funcdef.head in [:(=), :function]
    prototype, funcbody = funcdef.args

    if prototype.head == :call
        prototype, destructure = _destructure(prototype)
    elseif prototype.head == :where
        typeargs = prototype.args[2:end]
        prototype, destructure = _destructure(prototype.args[1])
        prototype = :($prototype where $((esc(t) for t in typeargs)...))
    end

    quote
        $prototype = begin
            $(destructure...)
            $(esc(funcbody))
        end
    end
end

# Utility type for handling several types of indexing operations

struct IndexMap
    data :: OrderedDict{Int, ArrayEvaluable{Int}}

    IndexMap(data::OrderedDict{Int, ArrayEvaluable{Int}}) = new(data)
    function IndexMap(indices::Pair...)
        new(OrderedDict{Int, ArrayEvaluable{Int}}(k => asarray(v) for (k,v) in indices))
    end
    function IndexMap(indices...)
        new(OrderedDict{Int, ArrayEvaluable{Int}}(k => asarray(v) for (k,v) in enumerate(indices) if !isa(v, Colon)))
    end
end

Base.get(self::IndexMap, i::Int, default) = get(self.data, i, default)
Base.getindex(map::IndexMap, array) = getindex(array, (get(map, i, :) for i in 1:ndims(array))...)
Base.getindex(self::IndexMap, i::Int) = get(self.data, i, :)
Base.isempty(self::IndexMap) = isempty(self.data)
Base.iterate(self::IndexMap) = iterate(self.data)
Base.iterate(self::IndexMap, state) = iterate(self.data, state)
Base.keys(self::IndexMap) = keys(self.data)
Base.ndims(self::IndexMap, i::Int) = i in keys(self.data) ? ndims(self.data[i]) : 1
Base.size(self::IndexMap, i::Int) = size(self.data[i])
Base.size(self::IndexMap, i::Int, j::Int) = size(self.data[i], j)
Base.size(self::IndexMap, s::Shape, i::Int) = i in keys(self.data) ? size(self.data[i]) : (s[i],)
Base.values(self::IndexMap) = values(self.data)

dimcheck(self::IndexMap, low::Int=0, hi::Int=typemax(Int)) =
    all(low <= ndims(self, i) <= hi for i in keys(self))

resultndims(self::IndexMap, d::Int) = reduce(+, (ndims(v)-1 for v in values(self.data)); init=0) + d
resultsize(self::IndexMap, shape::Shape) = Tuple(flatten(size(self, shape, i) for i in 1:length(shape)))

function codegen(self::IndexMap, values, d::Int)
    subst = OrderedDict(zip(keys(self), values))
    ((get(subst, i, :) for i in 1:d)...,)
end

# Miscellaneous

function insertmany!(target::Vector{T}, axes, value::T) where T
    sortaxes = reverse(sort(collect(Int, axes)))
    for ax in sortaxes
        insert!(target, ax, value)
    end
end

max_stars_repo_name: TheBB/Jutils

max_issues_repo_name: TheBB/Jutils

max_forks_repo_name: TheBB/Jutils


text:
using Revise
using ApproxFun, LinearAlgebra, Parameters, Setfield

using PseudoArcLengthContinuation, Plots
const PALC = PseudoArcLengthContinuation

####################################################################################################
# specific methods for ApproxFun
import Base: eltype, similar, copyto!, length
import LinearAlgebra: mul!, rmul!, axpy!, axpby!, dot, norm

similar(x::ApproxFun.Fun, T) = (copy(x))
similar(x::ApproxFun.Fun) = copy(x)
mul!(w::ApproxFun.Fun, v::ApproxFun.Fun, α) = (w .= α * v)

eltype(x::ApproxFun.Fun) = eltype(x.coefficients)
length(x::ApproxFun.Fun) = length(x.coefficients)

dot(x::ApproxFun.Fun, y::ApproxFun.Fun) = sum(x * y)

# do not put y .= a .* x .+ y, this puts a lot of coefficients!
axpy!(a::Float64, x::ApproxFun.Fun, y::ApproxFun.Fun) = (y .= a * x + y)
axpby!(a::Float64, x::ApproxFun.Fun, b::Float64, y::ApproxFun.Fun) = (y .= a * x + b * y)
rmul!(y::ApproxFun.Fun, b::Float64) = (y.coefficients .*= b; y)
rmul!(y::ApproxFun.Fun, b::Bool) = b == true ? y : (y.coefficients .*= 0; y)

# copyto!(x::ApproxFun.Fun, y::ApproxFun.Fun) = ( copyto!(x.coefficients, y.coefficients);x)
copyto!(x::ApproxFun.Fun, y::ApproxFun.Fun) = ( (x.coefficients = copy(y.coefficients);x))

####################################################################################################

N(x; a = 0.5, b = 0.01) = 1 + (x + a * x^2) / (1 + b * x^2)
dN(x; a = 0.5, b = 0.01) = (1 - b * x^2 + 2 * a * x)/(1 + b * x^2)^2

function F_chan(u, p)
	@unpack alpha, beta = p
	return [Fun(u(0.), domain(u)) - beta,
			Fun(u(1.), domain(u)) - beta,
			Δ * u + alpha * N(u, b = beta)]
end

function dF_chan(u, v, p)
	@unpack alpha, beta = p
	return [Fun(v(0.), domain(u)),
			Fun(v(1.), domain(u)),
			Δ * v + alpha * dN(u, b = beta) * v]
end

function Jac_chan(u, p)
	@unpack alpha, beta = p
	return [Evaluation(u.space, 0.),
			Evaluation(u.space, 1.),
			Δ + alpha * dN(u, b = beta)]
end

function finalise_solution(z, tau, step, contResult)
	printstyled(color=:red,"--> AF length = ", (z, tau) .|> length ,"\n")
	chop!(z.u, 1e-14);chop!(tau.u, 1e-14)
	true
end

sol = Fun( x -> x * (1-x), Interval(0.0, 1.0))
const Δ = Derivative(sol.space, 2);
par_af = (alpha = 3., beta = 0.01)

optnew = NewtonPar(tol = 1e-12, verbose = true)
	out, _, flag = @time PALC.newton(
		F_chan, Jac_chan, sol, par_af, optnew, normN = x -> norm(x, Inf64))
	# Plots.plot(out, label="Solution")

optcont = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, pMax = 4.1, plotEveryNsteps = 10, newtonOptions = NewtonPar(tol = 1e-8, maxIter = 20, verbose = true), maxSteps = 300)

	br, _ = @time continuation(
		F_chan, Jac_chan, out, par_af, (@lens _.alpha), optcont;
		plot = true,
		plotSolution = (x, p; kwargs...) -> plot!(x; label = "l = $(length(x))", kwargs...),
		verbosity = 2,
		normC = x -> norm(x, Inf64))
####################################################################################################
# Example with deflation technique
deflationOp = DeflationOperator(2.0, (x, y) -> dot(x, y), 1.0, [out])
par_def = @set par_af.alpha = 3.3

optdef = setproperties(optnew; tol = 1e-9, maxIter = 1000)

solp = copy(out)
	solp.coefficients .*= (1 .+ 0.41*rand(length(solp.coefficients)))

plot(out);plot!(solp)

outdef1, _, flag = @time newton(
	F_chan, Jac_chan,
	solp, par_def,
	optdef, deflationOp)
	flag && push!(deflationOp, outdef1)

plot(deflationOp.roots)
####################################################################################################
# other dot product
# dot(x::ApproxFun.Fun, y::ApproxFun.Fun) = sum(x * y) * length(x) # gives 0.1

optcont = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, pMax = 4.1, plotEveryNsteps = 10, newtonOptions = NewtonPar(tol = 1e-8, maxIter = 20, verbose = true), maxSteps = 300, theta = 0.2)

	br, _ = @time continuation(
		F_chan, Jac_chan, out, par_af, (@lens _.alpha), optcont;
		dotPALC = (x, y) -> dot(x, y),
		plot = true,
		# finaliseSolution = finalise_solution,
		plotSolution = (x, p; kwargs...) -> plot!(x; label = "l = $(length(x))", kwargs...),
		verbosity = 2,
		# printsolution = x -> norm(x, Inf64),
		normC = x -> norm(x, Inf64))
####################################################################################################
# tangent predictor with Bordered system
br, _ = @time continuation(
	F_chan, Jac_chan, out, par_af, (@lens _.alpha), optcont,
	tangentAlgo = BorderedPred(),
	plot = true,
	finaliseSolution = finalise_solution,
	plotSolution = (x, p;kwargs...)-> plot!(x; label = "l = $(length(x))", kwargs...))
####################################################################################################
# tangent predictor with Bordered system
# optcont = @set optcont.newtonOptions.verbose = true
indfold = 2
outfold, _, flag = @time newtonFold(
		F_chan, Jac_chan,
		br, indfold, #index of the fold point
		par_af, (@lens _.alpha),
		optcont.newtonOptions)
	flag && printstyled(color=:red, "--> We found a Fold Point at α = ", outfold[end], ", β = 0.01, from ", br.bifpoint[indfold][3],"\n")
#################################################################################################### Continuation of the Fold Point using minimally augmented
indfold = 2

outfold, _, flag = @time newtonFold(
	(x, p) -> F_chan(x, p),
	(x, p) -> Jac_chan(x, p),
	(x, p) -> Jac_chan(x, p),
	br, indfold, #index of the fold point
	optcont.newtonOptions)

ext: jl


text:
#Functions for the Harmonic correction in 2ls.
using LinearAlgebra
using QuadGK
using FastClosures

function correction_matrix(tf, env_Δ, env_pulse1, env_pulse2)
	nΔ = length(env_Δ)
	npulse1, npulse2 = length(env_pulse1), length(env_pulse2)
	dims = nΔ + npulse1 + npulse2
    M = zeros(3, dims)
    for i = 1: nΔ
        h1 = env_Δ[i]
        func1 = @closure t-> h1(t/tf)*sin(2*c(t, tf))
        int_func1, err = quadgk(func1, 0.0, tf, rtol=1e-10, atol=1e-14)
        func2 = @closure t-> h1(t/tf)*cos(2*c(t, tf))
        int_func2, err = quadgk(func2, 0.0, tf, rtol=1e-10, atol=1e-14)
        M[2, i], M[3, i] = int_func1, int_func2
    end
    for i = 1: npulse1
        h2 = env_pulse1[i]
        func3 = @closure t-> h2(t/tf)*(1 + cos(2*t))/2
        int_func3, err = quadgk(func3, 0.0, tf, rtol=1e-10, atol=1e-14)
        func4 = @closure t-> -h2(t/tf)*sin(2*t)*cos(2*c(t, tf))/2
        int_func4, err = quadgk(func4, 0.0, tf, rtol=1e-10, atol=1e-14)
        func5 = @closure t-> h2(t/tf)*sin(2*t)*sin(2*c(t, tf))/2
        int_func5, err = quadgk(func5, 0.0, tf, rtol=1e-10, atol=1e-14)

        col = nΔ + i
        M[1, col], M[2, col], M[3, col] = int_func3, int_func4, int_func5
    end
    for i = 1: npulse2
        h3 = env_pulse2[i]
        func6 = @closure t-> h3(t/tf)*sin(2*t)/2.0
        int_func6, err = quadgk(func6, 0.0, tf, rtol=1e-10, atol=1e-14)
        func7 = @closure t-> -h3(t/tf)*(1 - cos(2*t))*cos(2*c(t, tf))/2
        int_func7, err = quadgk(func7, 0.0, tf, rtol=1e-10, atol=1e-14)
        func8 = @closure t-> h3(t/tf)*(1 - cos(2*t))*sin(2*c(t, tf))/2
        int_func8, err = quadgk(func8, 0.0, tf, rtol=1e-10, atol=1e-14)

        col = nΔ + npulse1 + i
        M[1, col], M[2, col], M[3, col] = int_func6, int_func7, int_func8
    end
    return M
end

function coefficients_W(Omega, M)
    ih_vec = Omega[1: 3]
    coeffs = pinv(M)*ih_vec
    return coeffs
end

function g(t, tf, coeffs, env_pulse1, env_pulse2)
	npulse1, npulse2 = length(env_pulse1), length(env_pulse2)
	cost, sint = cos(t), sin(t)
	gt = zero(t)
    for i = 1: npulse1
        h2 = env_pulse1[i]
        gt += coeffs[i]*h2(t/tf)*cost
    end
    for i = 1: npulse2
        h3 = env_pulse2[i]
        col = npulse1 + i
        gt += coeffs[col]*h3(t/tf)*sint
    end
	return gt
end
function Δ(t, tf, coeffs, env_Δ)
	nΔ = length(env_Δ)
	Δt = zero(t)
    for i = 1: nΔ
        h1 = env_Δ[i]
        Δt += coeffs[i]*h1(t/tf)
    end
	return Δt
end
function Wharmonic(t, tf, coeffs, env_Δ, env_pulse1, env_pulse2)
	Δt = Δ(t, tf, coeffs[1], env_Δ)
	gt = g(t, tf, coeffs[2], env_pulse1, env_pulse2)

    W = zeros(4)
    W[1] = gt*cos(t)
    W[2] = -gt*sin(t)
    W[3] = Δt
    return W
end

function Wharmonic_I(t, tf, coeffs, env_Δ, env_pulse1, env_pulse2)
	cost, sint = cos(t), sin(t)
	cosθ, sinθ = cos(2*c(t, tf)), sin(2*c(t, tf))
	Δt = Δ(t, tf, coeffs[1], env_Δ)
	gt = g(t, tf, coeffs[2], env_pulse1, env_pulse2)

	W_I = zeros(4)
	W_I[1] = gt*cost
	W_I[2] = Δt*sinθ - gt*sint*cosθ
	W_I[3] = Δt*cosθ + gt*sint*sinθ
	return W_I
end

function get_coeffs(M, coeffs, order, tf, env_Δ, env_pulse1, env_pulse2)
	nΔ = length(env_Δ)
	A = @closure (t-> - V_I(t, tf)
					  - Wharmonic_I(t, tf, coeffs, env_Δ, env_pulse1, env_pulse2))

    Ndims = 4
    Omegaf = solveOmega(tf, order, A, comm!, Ndims)

    vec = sum(Omegaf, dims=2)
    coeffs = coefficients_W(vec, M)
    return (coeffs[1: nΔ], coeffs[nΔ + 1: end])
end
function get_coeffs(M, order, tf, env_Δ, env_pulse1, env_pulse2)
    nΔ = length(env_Δ)
	npulse1, npulse2 = length(env_pulse1), length(env_pulse2)
    coeffs = (zeros(nΔ),  zeros(npulse1 + npulse2))
    for k = 1: order
		coeffs2 = get_coeffs(M, coeffs, k, tf, env_Δ, env_pulse1, env_pulse2)
        coeffs = coeffs .+ coeffs2
    end
    return coeffs
end

function U(tf, σ, coeffs, env_Δ, env_pulse1, env_pulse2)
	H = @closure (t-> H0(t, tf, true)
					  + Wharmonic(t, tf, coeffs, env_Δ, env_pulse1, env_pulse2))

    sol = solveU(H, tf, σ)
    return sol(tf)
end

max_stars_repo_name: thalesfr/MagnusCorrection

max_issues_repo_name: thalesfr/MagnusCorrection

max_forks_repo_name: thalesfr/MagnusCorrection


text:
using UtilitiesForMRI, LinearAlgebra, CUDA, Test
CUDA.allowscalar(false)

# Cartesian domain
n = (256, 256, 256)
h = (abs(randn()), abs(randn()), abs(randn()))
X = spatial_sampling(Float64, n; h=h)

# Cartesian sampling in k-space
phase_encoding = (1,2)
K = kspace_Cartesian_sampling(X; phase_encoding=phase_encoding)

# Adjoint test (linear operator)
P = phase_shift(K)
nt, nk = size(K)
τ = randn(Float64, nt, 3)
Pτ = P(τ)
d = randn(ComplexF64, nt, nk)
e = randn(ComplexF64, nt, nk)
@test dot(Pτ*d, e) ≈ dot(d, Pτ'*e) rtol=1e-6

# Jacobian
Pτd, Pτ, ∂Pτd = ∂(P()*d, τ)
@test Pτd ≈ Pτ*d rtol=1e-6

# Adjoint test (Jacobian)
Δτ = randn(ComplexF64, nt, 3); Δτ *= norm(τ)/norm(Δτ)
Δd_ = ∂Pτd*Δτ; Δd = randn(ComplexF64, size(Δd_)); Δd *= norm(Δd_)/norm(Δd)
@test dot(∂Pτd*Δτ, Δd) ≈ dot(Δτ, ∂Pτd'*Δd) rtol=1e-6

# Gradient test
Δτ = randn(Float64, nt, 3); Δτ *= norm(τ)/norm(Δτ)
t = 1f-6
Pτd_p1 = P(τ+0.5*t*Δτ)*d
Pτd_m1 = P(τ-0.5*t*Δτ)*d
@test (Pτd_p1-Pτd_m1)/t ≈ ∂Pτd*Δτ rtol=1e-5

max_stars_repo_name: grizzuti/UtilitiesForMRI

max_issues_repo_name: grizzuti/UtilitiesForMRI

max_forks_repo_name: grizzuti/UtilitiesForMRI


text:
"""
    mat_center(X; mat = "her")

This function centers input.
If X is a Hermitian matrix, then ``Xcen = X - Tr(X)/N*I```.
If X is a rectangular matrix, then ``Xcen = X - X*ones(M,1)*ones(1,M)/M``.

# Arguments
-   `X`: a Hermitian or rectangular matrix
-   `mat`: the type of `X`, valid options are "her" and "rec"

# Outputs:
-   `Xc`: centered version of `X`
"""
function mat_center(X::Array{T,2}; mat = "her") where T <: Number
    N, M = size(X)

    if mat == "her"
        # for Hermitian matrix, X = X - I*tr(X)/N
        if M != N
            error("X should be a symmetric matrix")
        end
        Xc = X - I*tr(X)/M
    elseif mat == "rec"
        # for rectangular matrix, A = A - mean(A)*ones(N, M)
        Xc = X - mean(X)*ones(N, M)
    end

    return Xc
end

"""
    free_whiten(Z; mat = "her")

This function whitens input array of matrices in the free probability sense

# Arguments
-   `Z`: an array of matrices of mat and of the same dimensions
-   `mat`: the type of Z[i], valid options are "her" and "rec"

# Outputs
-   `Y`: an array of centered matrices, whitened in the free probability sense,
       that is, Tr(Y[i]*Y[j]')/size(Y[1],1) = (i == j)
-   `U`: a matrix of size s-by-s, where s = size(Z, 1)
-   `Σ`: a matrix of size s-by-s,
"""
function free_whiten(Z::Array{Array{T, 2}}; mat = "her") where T <: Number
    Zc = mat_center.(Z; mat = mat)
    N = size(Zc[1], 1)

    # compute the cov matrix such that
    # cov[i,j] = tr(Zc[i]*Zc[j]')/N
    eigcov = eigen(real(tr.(Zc*Zc'/N)));
    Λ, U = eigcov.values, eigcov.vectors # only take the real part

    Λ[Λ.<0] .= 0 # keep non-negative eigs
    Σ = Diagonal(Λ).^0.5
    Y = U*pinv(Σ)*U'*Zc

    return Y, U, Σ
end

# if Z is just a single matrix, we first convert it to [Y].
function free_whiten(Z::Array{T, 2}; mat = "her") where T <: Number
    return free_whiten([Z]; mat = mat)
end


text:
function newtonRaphson(f::Function, fprime::Function, xi::Number, ea::Number, maxIter::Int64)
  # Estimate roots of functions using the Newton-Raphson method
  for i = 1:maxIter
    xiprev  = xi;
    xi = xi - (f(xi) / fprime(xi));

    if abs(xi - xiprev) / xi < ea
      break;
    end;
  end
  return xi;
end

max_stars_repo_path: Julia/Root Estimation/newtonraphson.jl

max_stars_repo_name: alexjohnj/numerical-methods

max_issues_repo_path: Julia/Root Estimation/newtonraphson.jl

max_issues_repo_name: alexjohnj/numerical-methods

max_forks_repo_path: Julia/Root Estimation/newtonraphson.jl

max_forks_repo_name: alexjohnj/numerical-methods


text:
data = open("../../../externals/Core.Math.Samples/data/Pejcic_318.csv", "r") do f   # "r" for reading
    filecontent = read(f,String) # attention that it can be used only once. The second time, without reopening the file, read() would return an empty string
    print(filecontent)
  end


text:
using Mappings

using Test
v = collect(1:1:26)
v
D1 = Dict(zip('A':'Z',v))
D2 = Dict(zip(v,rand(26)))
m1 = AssociativeMap(D1)
m2 = AssociativeMap(D2)

import Base.HasLength
m1.amap
image(m1,'A')
codom(m1)
dom(m2)
dom(m2)
codom(m1)
inverse(m1)
inverse(inverse(m1))

m1('A')

m1 ∘ m2

image(d)

image(m1,'A')

inverse(m)
preimage(m1)

invert(m.amap)

D1
using BenchmarkTools
m.(Set(('A', 'B', 'C')))
Set(m.(['A', 'B', 'C']))
[(k, v) for (k, v) in collect(D1)]

dict = Dict(1:10 .=> repeat(1:5, 2))
rdict = Dict(values(dict) .=> keys(dict))

D = invD

invD

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


text:
const API_ROOT = "https://www.quandl.com/api/v3"
const SERIES_ROOT = string(API_ROOT, "/datasets")
const DATA_TABLE_ROOT = string(API_ROOT, "/datatables")

function enhance(url::AbstractString, key::AbstractString, value, check::Function)
    value == nothing && return url
    check(value) || throw(ArgumentError("Invalid value for key `$key`: $value"))
    return string(url, "&", key, "=", value)
end

positive_integer(x) = x isa Integer && x > 0

non_negative_integer(x) = x isa Integer && x >= 0

valid_date(x) = x isa Date ||
    (x isa AbstractString && match(r"\d\d\d\d-\d\d-\d\d", x) !== nothing)

valid_order_string(x) = x ∈ ["asc", "desc"]

valid_sample_frequency(x) = x ∈ ["none", "daily", "weekly", "monthly", "quarterly", "annual"]

valid_transform(x) = x ∈ ["none", "diff", "rdiff", "rdiff_from", "cumul", "normalize"]

function make_url(q::Quandl, series::TimeSeries;
             limit = nothing, column_index = nothing,
             start_date = nothing, end_date = nothing,
             order = nothing, collapse = nothing, transform = nothing)

    url = "$SERIES_ROOT/$(series.code).csv?api_key=$(q.api_key)"
    url = enhance(url, "limit", limit, positive_integer)
    url = enhance(url, "column_index", column_index, non_negative_integer)
    url = enhance(url, "start_date", start_date, valid_date)
    url = enhance(url, "end_date", end_date, valid_date)
    url = enhance(url, "order", order, valid_order_string)
    url = enhance(url, "collapse", collapse, valid_sample_frequency)
    url = enhance(url, "transform", transform, valid_transform)
    return url
end

# filter condition operators
eq(x, y) = FilterCondition("=", x, y)
gt(x, y) = FilterCondition(".gt=", x, y)
lt(x, y) = FilterCondition(".lt=", x, y)
gte(x, y) = FilterCondition(".gte=", x, y)
lte(x, y) = FilterCondition(".lte=", x, y)

function enhance(url::AbstractString, cond::FilterCondition)
    return string(url, "&", cond.field, cond.op, cond.value)
end

function make_url(q::Quandl, table::Table;
            filters::AbstractVector{T} = FilterCondition[],
            columns::AbstractVector{String} = String[]) where {T <: FilterCondition}

    url = "$DATA_TABLE_ROOT/$(table.code).csv?api_key=$(q.api_key)"
    for cond in filters
        url = enhance(url, cond)
    end
    if length(columns) > 0
        url = string(url, "&qopts.columns=", join(columns, ","))
    end
    return url
end

strip_key(url::AbstractString) = replace(url, r"api_key=[0-9a-zA-Z]*" => "api_key=xxxxxx")

function (q::Quandl)(x::T; verbose = false, kwargs...) where {T <: Union{TimeSeries, Table}}
    url = make_url(q, x; kwargs...)
    verbose && @info "Requesting data" strip_key(url)
    try
        response = HTTP.get(url)
        return DataFrame(CSV.File(response.body))
    catch ex
        @error "Unable to download data" ex
        return nothing
    end
end


text:
module grpcerrors
  const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
  include("grpcerrors_pb.jl")
end


text:
DEBUG = true
if DEBUG
    day = 9
else
    if length(ARGS) != 1
        quit()
    end
    day = strip(ARGS[1])
    day = parse( Int64, day)
end
g_year = "2017"
g_month = "09"

println("processing $(g_year)-$(g_month)-$(day)")
flush( STDOUT )
#################################### start processing
push!( LOAD_PATH, raw"/home/liye/project/code/lib_julia/vmim")
@time using vmim

### get whole_data_data list from json
f = open( @sprintf("/home/liye/working_paper/2018-forecasting/%s-%s-vmim-files.json", g_year, g_month) )
str = readstring( f )
files = JSON.parse( str );

## process file of one day
area_path = "/home/liye/project/code/2-near-missing/rect.json"
out_path = "/home/liye/project/data/2-near-missing/data/$(g_month)/" * string(day)
if !ispath(out_path)
    mkdir(out_path)
end
extract_in_rect_sampling( "2017", "$(g_month)", day, files, area_path, out_path )

max_stars_repo_name: juliatrend/website

max_issues_repo_name: juliatrend/website

max_forks_repo_name: juliatrend/website


text:
@inline function adjunij!(mat::AbstractMatrix, locs, U::Matrix)
    for j = 1:size(U, 2)
        for i = 1:size(U, 2)
            @inbounds U[i, j] += mat[locs[i], locs[j]]
        end
    end
    return U
end

@inline function adjunij!(mat::AbstractMatrix, locs, U::SparseMatrixCSC)
    @inbounds for j = 1:size(U, 2)
        S = U.colptr[j]
        E = U.colptr[j+1] - 1
        for ii = S:E
            @inbounds U.nzval[ii] += mat[locs[U.rowval[ii]], locs[j]]
        end
    end
    return U
end

@inline function adjunij!(mat::SDDiagonal, locs, U::Diagonal)
    @inbounds U.diag .+= mat.diag[locs]
    return U
end

@inline function adjunij!(mat::AbstractMatrix, locs, U::Diagonal)
    @inbounds for i = 1:size(U, 1)
        li = locs[i]
        U.diag[i] += mat[li, li]
    end
    return U
end

@inline function adjunij!(mat::SDPermMatrix, locs, U::PermMatrix)
    @inbounds U.vals .+= mat.vals[locs]
    return U
end

@inline function adjunij!(mat::AbstractMatrix, locs, U::PermMatrix)
    for i = 1:size(U, 1)
        @inbounds U.vals[i] += mat[locs[i], locs[U.perm[i]]]
    end
    return U
end

function adjcunmat(
    adjy::AbstractMatrix,
    nbit::Int,
    cbits::NTuple{C,Int},
    cvals::NTuple{C,Int},
    U0::AbstractMatrix{T},
    locs::NTuple{M,Int},
) where {C,M,T}
    U, ic, locs_raw = YaoBlocks.reorder_unitary(nbit, cbits, cvals, U0, locs)
    adjU = _render_adjU(U)

    controldo(ic) do i
        adjunij!(adjy, locs_raw .+ i, adjU)
    end

    adjU =
        all(diff([locs...]) .> 0) ? adjU :
        YaoBase.reorder(adjU, collect(locs) |> sortperm |> sortperm)
    adjU
end

@inline function adju1ij!(csc::SparseMatrixCSC{T}, i::Int, j::Int, adjU::Matrix) where {T}
    @inbounds begin
        adjU[1, 1] += csc.nzval[2*i-1]
        adjU[1, 2] += csc.nzval[2*i]
        adjU[2, 1] += csc.nzval[2*j-1]
        adjU[2, 2] += csc.nzval[2*j]
    end
    adjU
end

function adju1mat(adjy, nbit::Int, U1::SDMatrix, ibit::Int)
    mask = bmask(ibit)
    step = 1 << (ibit - 1)
    step_2 = 1 << ibit

    adjy = projection(YaoBlocks._initialize_output(nbit, 0, U1), adjy)
    adjU = _render_adjU(U1)

    for j = 0:step_2:1<<nbit-step
        @inbounds @simd for i = j+1:j+step
            adju1ij!(adjy, i, i + step, adjU)
        end
    end
    adjU
end

_render_adjU(U0::AbstractMatrix{T}) where {T} = zeros(T, size(U0)...)
_render_adjU(U0::SDSparseMatrixCSC{T}) where {T} = SparseMatrixCSC(
    size(U0)...,
    dynamicize(U0.colptr),
    dynamicize(U0.rowval),
    zeros(T, U0.nzval |> length),
)
_render_adjU(U0::SDDiagonal{T}) where {T} = Diagonal(zeros(T, size(U0, 1)))
_render_adjU(U0::SDPermMatrix{T}) where {T} = PermMatrix(U0.perm, zeros(T, length(U0.vals)))

# DEPRECATED
#=
@inline function adjsetcol!(csc::SparseMatrixCSC, icol::Int, rowval::AbstractVector, nzval::SubArray)
    begin
        S = csc.colptr[icol]
        E = csc.colptr[icol+1]-1
        nzval .+= view(csc.nzval, S:E)
    end
    csc
end

@inline function adjunij!(mat::SparseMatrixCSC, locs, U::Matrix)
    for j = 1:size(U, 2)
        adjsetcol!(mat, locs[j], locs, view(U,:,j))
    end
    return U
end

@inline function adjunij!(mat::SparseMatrixCSC, locs, U::SparseMatrixCSC)
    for j = 1:size(U, 2)
        S = U.colptr[j]
        E = U.colptr[j+1]-1
        @inbounds adjsetcol!(mat, locs[j], locs, view(U.nzval,S:E))
    end
    return U
end
=#


text:
using StatsFuns: logsumexp, log1pexp

using CUDA: CUDA, @cuda
using DataFrames: DataFrame
using LoopVectorization: @avx
using LogicCircuits: balance_threads, BitCircuit

export marginal, MAR, marginal_all, marginal_log_likelihood,
marginal_log_likelihood_avg, marginal_flows, marginal_flows_down

#####################
# Circuit marginal evaluation
#####################

"""
Evaluate a probabilistic circuit as a function
"""
function (root::ProbCircuit)(data...)
    marginal(root, data...)
end

"""
    marginal(root::ProbCircuit, data::Union{Real,Missing}...)
    marginal(root::ProbCircuit, data::Union{Vector{Union{Bool,Missing}},CuVector{UInt8}})
    marginal(circuit::ProbCircuit, data::DataFrame)
    marginal(circuit::ParamBitCircuit, data::DataFrame)::AbstractVector

Evaluate marginals of the circuit bottom-up for given input(s).

Missing values should be denoted by `missing` in the data.

Outputs ``\\log{p(x^o)}`` for each data point.
"""
marginal(root::ProbCircuit, data::Union{Real,Missing}...) =
    marginal(root, collect(Union{Bool,Missing}, data))

marginal(root::ProbCircuit, data::Union{Vector{Union{Bool,Missing}},CuVector{UInt8}}) =
    marginal(root, DataFrame(reshape(data, 1, :)))[1]

marginal(circuit::ProbCircuit, data::Union{DataFrame, Array{DataFrame}}) =
    marginal(same_device(ParamBitCircuit(circuit, data), data) , data)

function marginal(circuit::ParamBitCircuit, data::DataFrame)::AbstractVector
    marginal_all(circuit, data)[:,end]
end

function marginal(circuit::ParamBitCircuit, data::Array{DataFrame})::AbstractVector
    if isgpu(data)
        marginals = CuVector{Float64}(undef, num_examples(data))
    else
        marginals = Vector{Float64}(undef, num_examples(data))
    end

    v, start_idx = nothing, 1
    map(data) do d
        v = marginal_all(circuit, d, v)
        v_len = size(v, 1)
        @inbounds @views marginals[start_idx: start_idx + v_len - 1] .= v[:, end]
        start_idx += v_len
    end

    isgpu(data) ? convert(Array{Float64, 1}, to_cpu(marginals)) : marginals
end

"""
    MAR(pc, data)

Computes Marginal log likelhood of data. See docs for `marginal`.
"""
const MAR = marginal

"""
    marginal_log_likelihood(pc, data)

Compute the marginal likelihood of the PC given the data
"""
marginal_log_likelihood(pc, data) = begin
    if isweighted(data)
        # `data' is weighted according to its `weight' column
        data, weights = split_sample_weights(data)

        marginal_log_likelihood(pc, data, weights)
    else
        sum(marginal(pc, data))
    end
end
marginal_log_likelihood(pc, data, weights::DataFrame) = marginal_log_likelihood(pc, data, weights[:, 1])
marginal_log_likelihood(pc, data, weights::AbstractArray) = begin
    if isgpu(weights)
        weights = to_cpu(weights)
    end
    likelihoods = marginal(pc, data)
    if isgpu(likelihoods)
        likelihoods = to_cpu(likelihoods)
    end
    mapreduce(*, +, likelihoods, weights)
end
marginal_log_likelihood(pc, data::Array{DataFrame}; use_gpu::Bool = false) = begin
    if use_gpu
        data = to_gpu(data)
    end

    pbc = ParamBitCircuit(pc, data)

    total_ll::Float64 = 0.0
    if isweighted(data)
        data, weights = split_sample_weights(data)

        if use_gpu
            data = to_gpu(data)
        end

        v = nothing
        for idx = 1 : length(data)
            v = marginal_all(pbc, data[idx], v)
            likelihoods = v[:,end]
            if isgpu(likelihoods)
                likelihoods = to_cpu(likelihoods)
            end
            w = weights[idx]
            if isgpu(w)
                w = to_cpu(w)
            end
            total_ll += mapreduce(*, +, likelihoods, w)
        end
    else
        v = nothing
        for idx = 1 : length(data)
            v = marginal_all(pbc, data[idx], v)
            likelihoods = v[:,end]
            total_ll += sum(isgpu(likelihoods) ? to_cpu(likelihoods) : likelihoods)
        end
    end

    if use_gpu
        CUDA.unsafe_free!(v) # save the GC some effort
    end

    total_ll
end

"""
Compute the marginal likelihood of the PC given the data, averaged over all instances in the data
"""
marginal_log_likelihood_avg(pc, data) = begin
    if isweighted(data)
        # `data' is weighted according to its `weight' column
        data, weights = split_sample_weights(data)

        marginal_log_likelihood_avg(pc, data, weights)
    else
        marginal_log_likelihood(pc, data)/num_examples(data)
    end
end
marginal_log_likelihood_avg(pc, data, weights::DataFrame) = marginal_log_likelihood_avg(pc, data, weights[:, 1])
marginal_log_likelihood_avg(pc, data, weights) = begin
    if isgpu(weights)
        weights = to_cpu(weights)
    end
    marginal_log_likelihood(pc, data, weights)/sum(weights)
end
marginal_log_likelihood_avg(pc, data::Array{DataFrame}; use_gpu::Bool = false) = begin
    total_ll = marginal_log_likelihood(pc, data; use_gpu = use_gpu)

    if isweighted(data)
        data, weights = split_sample_weights(data)

        total_weights = 0.0
        for idx = 1 : length(weights)
            total_weights += sum(isgpu(weights[idx]) ? to_cpu(weights[idx]) : weights[idx])
        end

        avg_ll = total_ll / total_weights
    else
        avg_ll = total_ll / num_examples(data)
    end

    avg_ll
end

#####################
# Circuit evaluation of *all* nodes in circuit
#####################

"Evaluate the probabilistic circuit bottom-up for a given input and return the marginal probability value of all nodes"
marginal_all(circuit::ProbCircuit, data::DataFrame) =
    marginal_all(same_device(ParamBitCircuit(circuit, data), data) , data)

function marginal_all(circuit::ParamBitCircuit, data, reuse=nothing)
    @assert num_features(data) == num_features(circuit)
    values = init_marginal(data, reuse, num_nodes(circuit))
    marginal_layers(circuit, values)
    return values
end

"Initialize values from the data (data frames)"
function init_marginal(data, reuse, num_nodes; Float=Float32)
    flowtype = isgpu(data) ? CuMatrix{Float} : Matrix{Float}
    values = similar!(reuse, flowtype, num_examples(data), num_nodes)
    @views values[:,LogicCircuits.TRUE_BITS] .= log(one(Float))
    @views values[:,LogicCircuits.FALSE_BITS] .= log(zero(Float))
    # here we should use a custom CUDA kernel to extract Float marginals from bit vectors
    # for now the lazy solution is to move everything to the CPU and do the work there...
    data_cpu = to_cpu(data)
    nfeatures = num_features(data)
    for i=1:nfeatures
        marg_pos::Vector{Float} = log.(coalesce.(data_cpu[:,i], one(Float)))
        marg_neg::Vector{Float} = log.(coalesce.(1.0 .- data_cpu[:,i], one(Float)))
        values[:,2+i] .= same_device(marg_pos, values)
        values[:,2+nfeatures+i] .= same_device(marg_neg, values)
    end
    return values
end

# upward pass helpers on CPU

"Compute marginals on the CPU (SIMD & multi-threaded)"
function marginal_layers(circuit::ParamBitCircuit, values::Matrix)
    bc::BitCircuit = circuit.bitcircuit
    els = bc.elements
    pars = circuit.params
    for layer in bc.layers[2:end]
        Threads.@threads for dec_id in layer
            j = @inbounds bc.nodes[1,dec_id]
            els_end = @inbounds bc.nodes[2,dec_id]
            if j == els_end
                assign_marginal(values, dec_id, els[2,j], els[3,j], pars[j])
                j += 1
            else
                assign_marginal(values, dec_id, els[2,j], els[3,j], els[2,j+1], els[3,j+1], pars[j], pars[j+1])
                j += 2
            end
            while j <= els_end
                accum_marginal(values, dec_id, els[2,j], els[3,j], pars[j])
                j += 1
            end
        end
    end
end

assign_marginal(v::Matrix{<:AbstractFloat}, i, e1p, e1s, p1) =
    @views @. @avx v[:,i] = v[:,e1p] + v[:,e1s] + p1

accum_marginal(v::Matrix{<:AbstractFloat}, i, e1p, e1s, p1) = begin
    @avx for j=1:size(v,1)
        @inbounds x = v[j,i]
        @inbounds y = v[j,e1p] + v[j,e1s] + p1
        Δ = ifelse(x == y, zero(eltype(v)), abs(x - y))
        @inbounds  v[j,i] = max(x, y) + log1p(exp(-Δ))
    end
end

assign_marginal(v::Matrix{<:AbstractFloat}, i, e1p, e1s, e2p, e2s, p1, p2) = begin
    @avx for j=1:size(v,1)
        @inbounds x = v[j,e1p] + v[j,e1s] + p1
        @inbounds y = v[j,e2p] + v[j,e2s] + p2
        Δ = ifelse(x == y, zero(eltype(v)), abs(x - y))
        @inbounds  v[j,i] = max(x, y) + log1p(exp(-Δ))
    end
end

# upward pass helpers on GPU

"Compute marginals on the GPU"
function marginal_layers(circuit::ParamBitCircuit, values::CuMatrix;  dec_per_thread = 8, log2_threads_per_block = 8)
    circuit = to_gpu(circuit)
    bc = circuit.bitcircuit
    CUDA.@sync for layer in bc.layers[2:end]
        num_examples = size(values, 1)
        num_decision_sets = length(layer)/dec_per_thread
        num_threads =  balance_threads(num_examples, num_decision_sets, log2_threads_per_block)
        num_blocks = (ceil(Int, num_examples/num_threads[1]),
                      ceil(Int, num_decision_sets/num_threads[2]))
        @cuda threads=num_threads blocks=num_blocks marginal_layers_cuda(layer, bc.nodes, bc.elements, circuit.params, values)
    end
end

"CUDA kernel for circuit evaluation"
function marginal_layers_cuda(layer, nodes, elements, params, values)
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y
    for j = index_x:stride_x:size(values,1)
        for i = index_y:stride_y:length(layer)
            decision_id = @inbounds layer[i]
            k = @inbounds nodes[1,decision_id]
            els_end = @inbounds nodes[2,decision_id]
            @inbounds x = values[j, elements[2,k]] + values[j, elements[3,k]] + params[k]
            while k < els_end
                k += 1
                @inbounds y = values[j, elements[2,k]] + values[j, elements[3,k]] + params[k]
                Δ = ifelse(x == y, zero(eltype(values)), CUDA.abs(x - y))
                x = max(x, y) + CUDA.log1p(CUDA.exp(-Δ))
            end
            values[j, decision_id] = x
        end
    end
    return nothing
end

#####################
# Bit circuit marginals and flows (up and downward pass)
#####################

"Compute the marginal and flow of each node"
function marginal_flows(circuit::ProbCircuit, data,
    reuse_values=nothing, reuse_flows=nothing; on_node=noop, on_edge=noop, weights=nothing)
    bc = same_device(ParamBitCircuit(circuit, data), data)
    marginal_flows(bc, data, reuse_values, reuse_flows; on_node, on_edge, weights)
end

function marginal_flows(circuit::ParamBitCircuit, data,
            reuse_values=nothing, reuse_flows=nothing; on_node=noop, on_edge=noop, weights=nothing)
    @assert isgpu(data) == isgpu(circuit) "ParamBitCircuit and data need to be on the same device"
    values = marginal_all(circuit, data, reuse_values)
    flows = marginal_flows_down(circuit, values, reuse_flows; on_node, on_edge, weights)
    return values, flows
end

#####################
# Bit circuit marginal flows downward pass
#####################

"When marginals of nodes have already been computed, do a downward pass computing the marginal flows at each node"
function marginal_flows_down(circuit::ParamBitCircuit, values, reuse=nothing; on_node=noop, on_edge=noop, weights=nothing)
    flows = similar!(reuse, typeof(values), size(values)...)
    marginal_flows_down_layers(circuit, flows, values, on_node, on_edge; weights = weights)
    return flows
end

# downward pass helpers on CPU

"Evaluate marginals of the layers of a parameter bit circuit on the CPU (SIMD & multi-threaded)"
function marginal_flows_down_layers(pbc::ParamBitCircuit, flows::Matrix, values::Matrix, on_node, on_edge; weights = nothing)
    @assert flows !== values
    circuit::BitCircuit = pbc.bitcircuit
    els = circuit.elements
    for layer in Iterators.reverse(circuit.layers)
        Threads.@threads for dec_id in layer
            par_start = @inbounds circuit.nodes[3,dec_id]
            if iszero(par_start)
                if dec_id == num_nodes(circuit)
                    # marginal flow start from 0.0
                    @inbounds @views flows[:, dec_id] .= zero(eltype(flows))
                end
                # no parents, ignore (can happen for false/true node and root)
            else
                par_end = @inbounds circuit.nodes[4,dec_id]
                for j = par_start:par_end
                    par = @inbounds circuit.parents[j]
                    grandpa = @inbounds els[1,par]
                    sib_id = sibling(els, par, dec_id)
                    single_child = has_single_child(circuit.nodes, grandpa)
                    if single_child
                        if j == par_start
                            @inbounds @views @. flows[:, dec_id] = flows[:, grandpa]
                        else
                            accum_marg_flow(flows, dec_id, grandpa)
                        end
                    else
                        θ = eltype(flows)(pbc.params[par])
                        if j == par_start
                            assign_marg_flow(flows, values, dec_id, grandpa, sib_id, θ)
                        else
                            accum_marg_flow(flows, values, dec_id, grandpa, sib_id, θ)
                        end
                    end
                    # report edge flow only once:
                    sib_id > dec_id && on_edge(flows, values, dec_id, sib_id, par, grandpa, single_child, weights)
                end
            end
            on_node(flows, values, dec_id, weights)
        end
    end
end

function assign_marg_flow(f::Matrix{<:AbstractFloat}, v, d, g, s, θ)
    @inbounds @simd for j=1:size(f,1) #@avx gives incorrect results
        edge_flow = v[j, s] + v[j, d] - v[j, g] + f[j, g] + θ
        edge_flow = ifelse(isnan(edge_flow), typemin(eltype(f)), edge_flow)
        f[j, d] = edge_flow
    end
    # @assert !any(isnan, f[:,d])
end

function accum_marg_flow(f::Matrix{<:AbstractFloat}, d, g)
    @inbounds @simd for j=1:size(f,1) #@avx gives incorrect results
        x = f[j, d]
        y = f[j, g]
        Δ = ifelse(x == y, zero(eltype(f)), abs(x - y))
        f[j, d] = max(x, y) + log1p(exp(-Δ))
    end
    # @assert !any(isnan, f[:,d])
end

function accum_marg_flow(f::Matrix{<:AbstractFloat}, v, d, g, s, θ)
    @inbounds @simd for j=1:size(f,1) #@avx gives incorrect results
        x = f[j, d]
        y = v[j, s] + v[j, d] - v[j, g] + f[j, g] + θ
        y = ifelse(isnan(y), typemin(eltype(f)), y)
        Δ = ifelse(x == y, zero(eltype(f)), abs(x - y))
        f[j, d] = max(x, y) + log1p(exp(-Δ))
    end
    # @assert !any(isnan, f[:,d])
end

# downward pass helpers on GPU

"Pass marginal flows down the layers of a bit circuit on the GPU"
function marginal_flows_down_layers(pbc::ParamBitCircuit, flows::CuMatrix, values::CuMatrix,
            on_node, on_edge;
            dec_per_thread = 8, log2_threads_per_block = 7, weights = nothing)
    bc = pbc.bitcircuit
    CUDA.@sync for layer in Iterators.reverse(bc.layers)
        num_examples = size(values, 1)
        num_decision_sets = length(layer)/dec_per_thread
        num_threads =  balance_threads(num_examples, num_decision_sets, log2_threads_per_block)
        num_blocks = (ceil(Int, num_examples/num_threads[1]),
                      ceil(Int, num_decision_sets/num_threads[2]))
        @cuda threads=num_threads blocks=num_blocks marginal_flows_down_layers_cuda(layer, bc.nodes, bc.elements, bc.parents, pbc.params, flows, values, on_node, on_edge, weights)
    end
end

"CUDA kernel for passing marginal flows down circuit"
function marginal_flows_down_layers_cuda(layer, nodes, elements, parents, params, flows, values, on_node, on_edge, weights::Nothing)
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y
    for k = index_x:stride_x:size(values,1)
        for i = index_y:stride_y:length(layer)
            dec_id = @inbounds layer[i]
            if dec_id == size(nodes,2)
                # populate root flows
                flow = zero(eltype(flows))
            else
                par_start = @inbounds nodes[3,dec_id]
                flow = typemin(eltype(flows)) # log(0)
                if !iszero(par_start)
                    par_end = @inbounds nodes[4,dec_id]
                    for j = par_start:par_end
                        par = @inbounds parents[j]
                        grandpa = @inbounds elements[1,par]
                        v_gp = @inbounds values[k, grandpa]
                        prime = elements[2,par]
                        sub = elements[3,par]
                        θ = eltype(flows)(params[par])
                        if !iszero(v_gp) # edge flow only gets reported when non-zero
                            f_gp = @inbounds flows[k, grandpa]
                            single_child = has_single_child(nodes, grandpa)
                            if single_child
                                edge_flow = f_gp
                            else
                                v_prime = @inbounds values[k, prime]
                                v_sub = @inbounds values[k, sub]
                                edge_flow = compute_marg_edge_flow(v_prime, v_sub, v_gp, f_gp, θ)
                            end
                            flow = logsumexp_cuda(flow, edge_flow)
                            # report edge flow only once:
                            dec_id == prime && on_edge(flows, values, prime, sub, par, grandpa, k, edge_flow, single_child, nothing)
                        end
                    end
                end
            end
            @inbounds flows[k, dec_id] = flow
            on_node(flows, values, dec_id, k, flow, nothing)
        end
    end
    return nothing
end
function marginal_flows_down_layers_cuda(layer, nodes, elements, parents, params, flows, values, on_node, on_edge, weights::Any)
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y
    for k = index_x:stride_x:size(values,1)
        weight = @inbounds weights[k]
        for i = index_y:stride_y:length(layer)
            dec_id = @inbounds layer[i]
            if dec_id == size(nodes,2)
                # populate root flows
                flow = zero(eltype(flows))
            else
                par_start = @inbounds nodes[3,dec_id]
                flow = typemin(eltype(flows)) # log(0)
                if !iszero(par_start)
                    par_end = @inbounds nodes[4,dec_id]
                    for j = par_start:par_end
                        par = @inbounds parents[j]
                        grandpa = @inbounds elements[1,par]
                        v_gp = @inbounds values[k, grandpa]
                        prime = elements[2,par]
                        sub = elements[3,par]
                        θ = eltype(flows)(params[par])
                        if !iszero(v_gp) # edge flow only gets reported when non-zero
                            f_gp = @inbounds flows[k, grandpa]
                            single_child = has_single_child(nodes, grandpa)
                            if single_child
                                edge_flow = f_gp
                            else
                                v_prime = @inbounds values[k, prime]
                                v_sub = @inbounds values[k, sub]
                                edge_flow = compute_marg_edge_flow(v_prime, v_sub, v_gp, f_gp, θ, weight)
                            end
                            flow = logsumexp_cuda(flow, edge_flow)
                            # report edge flow only once:
                            dec_id == prime && on_edge(flows, values, prime, sub, par, grandpa, k, edge_flow, single_child, weight)
                        end
                    end
                end
            end
            @inbounds flows[k, dec_id] = flow
            on_node(flows, values, dec_id, k, flow, weight)
        end
    end
    return nothing
end

@inline function compute_marg_edge_flow(p_up, s_up, n_up, n_down, θ)
    x = p_up + s_up - n_up + n_down + θ
    ifelse(isnan(x), typemin(n_down), x)
end
@inline function compute_marg_edge_flow(p_up, s_up, n_up, n_down, θ, weight)
    x = p_up + s_up - n_up + n_down + θ
    ifelse(isnan(x), typemin(n_down), x)
end


text:
using Samsara, Conception
using Printf

function run(N::Int; number_of_nodes = nothing)
    case = init(number_of_nodes)
    demo(N, case)

    println("\n\n\n")
    all_data = []
    for it in case._all_nodes
        println("____________"*string(it)*"_____________________________________________")
        show_links_to(it, case)
        push!(all_data, get_maxvalue_histogram_to(case._all_nodes, it))
    end
    # all_data kan no brukes vidare, f.eks. ved å kalle `julia -i`
    all_data
end

function demo(N::Int, case)
    @assert N>0
    all_actions = case._all_actions
    all_nodes = case._all_nodes
    for i in 1:N
        _simulation_step!(case, all_actions)
    end
end

function _simulation_step!(case, all_actions)
    next_action = rand(all_actions)
    Samsara.set_action_in!(case, next_action)
    step_system_mechanics!(case)
end

function init(theN =nothing)
    muex = MuExS()
    all_actions=(SAT(:_up__), SAT(:down_), SAT(:noop_)) #, SAT(:reset))
    all_nodes = SAT[]
    if isnothing(theN)
        all_nodes  = [ SAT(:a, inMuExS=muex),
                       SAT(:b, inMuExS=muex),
                       SAT(:c, inMuExS=muex),
                       SAT(:d, inMuExS=muex),
                       SAT(:e, inMuExS=muex),
                       SAT(:f, inMuExS=muex),
                       SAT(:g, inMuExS=muex),
                       SAT(:h, inMuExS=muex),
                       SAT(:i, inMuExS=muex),
                     ]
    elseif isa(theN, Int)
        for i in 1:theN
            theName = @sprintf("node_%2d",i-1)
            theSAT = SAT(theName, inMuExS=muex)
            push!(all_nodes, theSAT)
        end
        for it in all_nodes
            @show it
        end
    else
        throw(ArgumentError("Wrong argument to init(theN)"))
    end
    case = Samsara.CircularSys(nodes= all_nodes, actions=all_actions )
end

function show_links_to(to_node::SAT, in_environment::Samsara.AbstractSystem)
    all_nodes = in_environment._all_nodes

    lista = sort(to_node._incoming_inferences, rev=true)
    for node in all_nodes
        println("__________________"*string(node)*" ⟶ "*string(to_node)*"__________________")
        show_all_elem_of(lista, node)
    end
end
function show_all_elem_of(lista::Vector, element)
    for it in lista
        (it.nodeA == element) && @show it
    end
end

# lag en maksvektor:
function get_maxvalue_histogram_to(node_lista::Vector, to_node::SAT)
    histrogram = Union{Missing, Float64}[]
    for it in node_lista
        push!(histrogram, get_maxvalue_from(to_node._incoming_inferences, it))
    end
    histrogram
end
function get_maxvalue_from(lista::Vector, node::SAT)
    all_elements = []
    for it in lista
        (it.nodeA == node) && push!(all_elements, inference(it))
    end
    sort!(all_elements)

    if length(all_elements) > 0
    @show all_elements
        return all_elements[end]
    else
        return missing
    end
end

println("\n************************ 200 *************************")
data1K = run(1000)  # ... med manuellt skrevet noder: a:j
println("\n************************ 1000 ****************************")
data1K = run(1000, number_of_nodes=5)
println("\n************************ 15000 ****************************")
data20K = run(15000, number_of_nodes=100)
#println("\n************************ 500.000 *************************")
#@show data100K = run(100000)

max_stars_repo_name: leikanger/Samsara

max_issues_repo_name: leikanger/Samsara

max_forks_repo_name: leikanger/Samsara


text:
using SnakePlot
using Test

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


text:
include("source.jl")

function init(S::Spins)
    S.N = para.N
    S.M = para.M

    S.A = rand(0,2, (para.N+2, para.M+2))
    S.A[1,:] = 0
	  S.A[para.N + 1, :] = 0
    S.A[:,1] = 0
	  S.A[:,para.M + 1] = 0

    boolarr1 = S.A[2:para.N + 2,2:para.M + 2] == 0
    rand = rand(para.N, para.M)
    boolarr2 = rand > 0.7
    boolarr = boolarr1 & boolarr2
    S.A[2:para.N + 2,2:para.M + 2] -= boolarr
    boolarr = boolarr1 & ~boolarr2
    S.A[2:para.N + 2,2:para.M + 2] += boolarr
end

max_stars_repo_name: anirbrhm/curious-eavesdropper

max_issues_repo_name: anirbrhm/curious-eavesdropper

max_forks_repo_name: anirbrhm/curious-eavesdropper


text:
export Square

"""
    Square{B} <: BCLattice{B, 2}

General square lattice with boundary condition `B`.

    Square([Boundary], height, width)
    Square([Boundary], shape)

Construct a square lattice. The shape follows the order of Julia's
Native Array.
"""
struct Square{B <: Boundary} <: BCLattice{B, 2}
    height::Int
    width::Int
end

Square(height::Integer, width::Integer) = Square{Fixed}(height, width)
Square(::Type{B}, height::Integer, width::Integer) where {B <: Boundary} =
    Square{B}(height, width)
Square(shape::Tuple{Int, Int}) = Square(Fixed, shape)
Square(::Type{B}, shape::Tuple{Int, Int}) where {B <: Boundary} =
    Square{B}(shape...)

shape(square::Square) = (square.height, square.width)
length(square::Square) = square.height * square.width
sites(square::Square) = SquareSiteIter(square)

struct SquareSiteIter{B<:Boundary} <: SiteIterator{Square{B}}
    height::Int
    width::Int
end

SquareSiteIter(square::Square{B}) where B = SquareSiteIter{B}(square.height, square.width)

struct SquareBondIter{B<:Boundary, K} <: BondIterator{Square{B}}
    height::Int
    width::Int
end

SquareBondIter(square::Square{B}, ::Type{K}) where {B, K} =
    SquareBondIter{B, K}(square.height, square.width)

lattice(itr::LatticeIterator{Square{B}}) where B =
    Square(B, itr.height, itr.width)

# TODO: add bounds check, throw LoadError when bond is invalid

function bonds(square::Square, bond::Integer)
    if bond % 2 == 1
        k = Int((bond + 1) / 2)
        return fuse(
            SquareBondIter(square, Vertical{k}),
            SquareBondIter(square, Horizontal{k})
        )
    else
        k = Int(bond / 2)
        return fuse(
            SquareBondIter(square, UpRight{k}),
            SquareBondIter(square, UpLeft{k})
        )
    end
end
#####################
# Iterator Interface
#####################

import Base: start, next, done, eltype, length, size, iteratorsize, HasShape

start(itr::SquareSiteIter) = (1, 1, 1)

function next(itr::SquareSiteIter, state)
    i, j, count = state

    if i > itr.height
        return (1, j+1), (2, j+1, count+1)
    end
    return (i, j), (i+1, j, count+1)
end

done(itr::SquareSiteIter, state) = state[3] > length(itr)

length(itr::SquareSiteIter) = itr.width * itr.height
size(itr::SquareSiteIter) = (itr.height, itr.width)
eltype(itr::SquareSiteIter) = Tuple{Int, Int}
iteratorsize(itr::SquareSiteIter) = HasShape()

######################
# Bonds
######################

start(itr::SquareBondIter) = (1, 1, 1)
eltype(itr::SquareBondIter) = NTuple{2, NTuple{2, Int}}
iteratorsize(itr::SquareBondIter) = HasShape()

function done(itr::SquareBondIter, state::NTuple{3, Int})
    i, j, count = state
    return count > length(itr)
end

##################
# Fixed Boundary
##################

length(itr::SquareBondIter{Fixed, Vertical{K}}) where K =
    (itr.height - K) * itr.width
length(itr::SquareBondIter{Fixed, Horizontal{K}}) where K =
    (itr.width - K) * itr.height
length(itr::SquareBondIter{Fixed, UpRight{K}}) where K =
    (itr.width - K) * (itr.height - K)
length(itr::SquareBondIter{Fixed, UpLeft{K}}) where K =
    (itr.width - K) * (itr.height - K)

size(itr::SquareBondIter{Fixed, Vertical{K}}) where K = (itr.height - K, itr.width)
size(itr::SquareBondIter{Fixed, Horizontal{K}}) where K = (itr.height, itr.width - K)
size(itr::SquareBondIter{Fixed, UpRight{K}}) where K = (itr.height - K, itr.width - K)
size(itr::SquareBondIter{Fixed, UpLeft{K}}) where K = (itr.height - K, itr.width - K)

function next(itr::SquareBondIter{Fixed, Vertical{K}}, state) where K
    i, j, count = state
    if i + K > itr.height
        return ((1, j+1), (1+K, j+1)), (2, j+1, count+1)
    end

    return ((i, j), (i+K, j)), (i+1, j, count+1)
end

function next(itr::SquareBondIter{Fixed, Horizontal{K}}, state) where K
    i, j, count = state
    if i > itr.height
        return ((1, j+1), (1, j+K+1)), (2, j+1, count+1)
    end
    return ((i, j), (i, j+K)), (i+1, j, count+1)
end

function next(itr::SquareBondIter{Fixed, UpRight{K}}, state) where K
    i, j, count = state
    if i+K > itr.height
        return ((1, j+1), (1+K, j+K+1)), (2, j+1, count+1)
    end
    return ((i, j), (i+K, j+K)), (i+1, j, count+1)
end

function next(itr::SquareBondIter{Fixed, UpLeft{K}}, state) where K
    i, j, count = state
    if i+K > itr.height
        return ((1+K, j+1), (1, j+K+1)), (2, j+1, count+1)
    end
    return ((i+K, j), (i, j+K)), (i+1, j, count+1)
end

####################
# Periodic Boundary
####################

length(itr::SquareBondIter{Periodic}) = itr.height * itr.width
size(itr::SquareBondIter{Periodic}) = (itr.height, itr.width)

function next(itr::SquareBondIter{Periodic, Vertical{K}}, state) where K
    i, j, count = state
    if i > itr.height
        return ((1, j+1), (K%itr.height+1, j+1)), (2, j+1, count+1)
    end

    return ((i, j), ((i+K-1)%itr.height+1, j)), (i+1, j, count+1)
end

function next(itr::SquareBondIter{Periodic, Horizontal{K}}, state) where K
    i, j, count = state
    if i > itr.height
        return ((1, j+1), (1, (j+K)%itr.width+1)), (2, j+1, count+1)
    end
    return ((i, j), (i, (j+K-1)%itr.width+1)), (i+1, j, count+1)
end

function next(itr::SquareBondIter{Periodic, UpRight{K}}, state) where K
    i, j, count = state
    if i > itr.height
        return ((1, j+1), (K%itr.height+1, (j+K)%itr.width+1)), (2, j+1, count+1)
    end
    return ((i, j), ((i+K-1)%itr.height+1, (j+K-1)%itr.width+1)), (i+1, j, count+1)
end

function next(itr::SquareBondIter{Periodic, UpLeft{K}}, state) where K
    i, j, count = state
    if i > itr.height
        return ((K%itr.height+1, j+1), (1, (j+K)%itr.width+1)), (2, j+1, count+1)
    end
    return (((i+K-1)%itr.height+1, j), (i, (j+K-1)%itr.width+1)), (i+1, j, count+1)
end

getid(square::Square, x::Int, y::Int) = x + (y - 1) * square.height


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    fit(estimator, X, z)

Fit estimator to coordinates `X` and values `z`,
and return a fitted estimator.
"""
function fit end

"""
    predict(estimator, xₒ)

Predict with the `estimator` at coordinates `xₒ`.
"""
function predict end

"""
    status(fitted)

Return the status of the `fitted` estimator. (e.g. the
factorization of the linear system was successful)
"""
function status end


text:
#=
    This code was generated by heta-compiler 0.6.6

=#

__platform__ = (function()

### create default constants
nameless_constants_num_ = NamedTuple{(

)}(Float64[

])

### create default observables
nameless_records_output_ = NamedTuple{(

)}(Bool[

])

### create default events
nameless_events_active_ = NamedTuple{(

)}(Bool[

])

### initialization of ODE variables and Records
function nameless_init_func_(cons)
    #() = cons

    # Heta initialize
    t = 0.0 # initial time

    # save results

    return ([ 0.0 ], # init for fake variable
        Float64[
        ]
    )
end

### calculate RHS of ODE
function nameless_ode_func_(du, u, p, t)
    cons = p.constants
    () = p.static
    () = u

    # Heta rules

    #p.static .= []
    du .= [ 1.0 ] # RHS for fake variable
end

### output function
function nameless_saving_generator_(outputIds::Vector{Symbol})
    function saving_(u, t, integrator)
        cons = integrator.p.constants
        () = integrator.p.static
        () = u

        # Heta rules

        # force amount

        d = Base.@locals
        return [d[id] for id in outputIds]
    end
end

### TIME EVENTS ###

### D EVENTS ###

### STOP EVENTS ###

### event assignments

### MODELS ###

nameless_model_ = (
  nameless_init_func_,
  nameless_ode_func_,
  NamedTuple{(
  )}([
  ]),
  NamedTuple{(
  )}([
  ]),
  NamedTuple{(
  )}([
  ]),
  nameless_saving_generator_,
  nameless_constants_num_,
  nameless_events_active_,
  nameless_records_output_
)
### create default constants
mod1_constants_num_ = NamedTuple{(
  :k_a,:Vd,:Vmax,:Km,:body_weight,:dose_first,:dose_last,:sigma_add,:sigma_prop,
)}(Float64[
  0.062,0.457,0.136,0.096,70,0.186,0.112,0.1,0,
])

### create default observables
mod1_records_output_ = NamedTuple{(
  :gut,:blood,:q_a,:C,:BrAC,:vabs,:vel,:sigma1,
)}(Bool[
  false,false,true,true,true,false,false,true,
])

### create default events
mod1_events_active_ = NamedTuple{(
  :adm1,:adm2,:adm3,:adm4,
)}(Bool[
  true,true,true,true,
])

### initialization of ODE variables and Records
function mod1_init_func_(cons)
    #(k_a,Vd,Vmax,Km,body_weight,dose_first,dose_last,sigma_add,sigma_prop,) = cons

    # Heta initialize
    t = 0.0 # initial time
    gut = 1.0
    C = 0.0
    BrAC = C / 1.0
    q_a = 0.0
    vabs = cons[1] * q_a
    blood = cons[2] * cons[5]
    vel = cons[3] * C / (cons[4] + C) * blood
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    # save results

    return (
        Float64[
            q_a,
            C * blood,
        ],
        Float64[
            gut,
            blood,
        ]
    )
end

### calculate RHS of ODE
function mod1_ode_func_(du, u, p, t)
    cons = p.constants
    (gut,blood,) = p.static
    (q_a,C_,) = u

    # Heta rules
    C = C_ / blood
    BrAC = C / 1.0
    vabs = cons[1] * q_a
    vel = cons[3] * C / (cons[4] + C) * blood
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    #p.static .= [gut,blood,]
    du .= [
      -vabs,  # dq_a_/dt
      vabs-vel,  # dC_/dt
    ]
end

### output function
function mod1_saving_generator_(outputIds::Vector{Symbol})
    function saving_(u, t, integrator)
        cons = integrator.p.constants
        (gut,blood,) = integrator.p.static
        (q_a,C_,) = u

        # Heta rules
        C = C_ / blood
        BrAC = C / 1.0
        vabs = cons[1] * q_a
        vel = cons[3] * C / (cons[4] + C) * blood
        sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

        # force amount
        q_a_ = q_a

        d = Base.@locals
        return [d[id] for id in outputIds]
    end
end

### TIME EVENTS ###
function mod1_adm1_tstops_func_(cons, times)
    return [0.0]
end
function mod1_adm2_tstops_func_(cons, times)
    return [0.333]
end
function mod1_adm3_tstops_func_(cons, times)
    return [0.667]
end
function mod1_adm4_tstops_func_(cons, times)
    return [1.5]
end

### D EVENTS ###

### STOP EVENTS ###

### event assignments
function mod1_adm1_affect_func_(integrator)
    cons = integrator.p.constants
    t = integrator.t
    (gut,blood,) = integrator.p.static
    (q_a,C_,) = integrator.u

    # Heta rules
    C = C_ / blood
    BrAC = C / 1.0
    vabs = cons[1] * q_a
    vel = cons[3] * C / (cons[4] + C) * blood
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    # recalculated values
    (q_a,) = (q_a + cons[6] * cons[5],)

    # save results
    integrator.u[[1,]] = [
        q_a,
    ]
    integrator.p.static[[]] = [
    ]
end
function mod1_adm2_affect_func_(integrator)
    cons = integrator.p.constants
    t = integrator.t
    (gut,blood,) = integrator.p.static
    (q_a,C_,) = integrator.u

    # Heta rules
    C = C_ / blood
    BrAC = C / 1.0
    vabs = cons[1] * q_a
    vel = cons[3] * C / (cons[4] + C) * blood
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    # recalculated values
    (q_a,) = (q_a + cons[6] * cons[5],)

    # save results
    integrator.u[[1,]] = [
        q_a,
    ]
    integrator.p.static[[]] = [
    ]
end
function mod1_adm3_affect_func_(integrator)
    cons = integrator.p.constants
    t = integrator.t
    (gut,blood,) = integrator.p.static
    (q_a,C_,) = integrator.u

    # Heta rules
    C = C_ / blood
    BrAC = C / 1.0
    vabs = cons[1] * q_a
    vel = cons[3] * C / (cons[4] + C) * blood
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    # recalculated values
    (q_a,) = (q_a + cons[6] * cons[5],)

    # save results
    integrator.u[[1,]] = [
        q_a,
    ]
    integrator.p.static[[]] = [
    ]
end
function mod1_adm4_affect_func_(integrator)
    cons = integrator.p.constants
    t = integrator.t
    (gut,blood,) = integrator.p.static
    (q_a,C_,) = integrator.u

    # Heta rules
    C = C_ / blood
    BrAC = C / 1.0
    vabs = cons[1] * q_a
    vel = cons[3] * C / (cons[4] + C) * blood
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    # recalculated values
    (q_a,) = (q_a + cons[7] * cons[5],)

    # save results
    integrator.u[[1,]] = [
        q_a,
    ]
    integrator.p.static[[]] = [
    ]
end

### MODELS ###

mod1_model_ = (
  mod1_init_func_,
  mod1_ode_func_,
  NamedTuple{(
    :adm1,
    :adm2,
    :adm3,
    :adm4,
  )}([
    (mod1_adm1_tstops_func_, mod1_adm1_affect_func_, :adm1, false),
    (mod1_adm2_tstops_func_, mod1_adm2_affect_func_, :adm2, false),
    (mod1_adm3_tstops_func_, mod1_adm3_affect_func_, :adm3, false),
    (mod1_adm4_tstops_func_, mod1_adm4_affect_func_, :adm4, false),
  ]),
  NamedTuple{(
  )}([
  ]),
  NamedTuple{(
  )}([
  ]),
  mod1_saving_generator_,
  mod1_constants_num_,
  mod1_events_active_,
  mod1_records_output_
)
### create default constants
mod2_constants_num_ = NamedTuple{(
  :k_a,:Vd,:Vmax,:Km,:body_weight,:dose_first,:dose_last,:sigma_add,:sigma_prop,:ke,
)}(Float64[
  0.062,0.457,0.136,0.096,70,0.186,0.112,0.1,0,0.011,
])

### create default observables
mod2_records_output_ = NamedTuple{(
  :gut,:blood,:q_a,:C,:BrAC,:vabs,:vel,:sigma1,
)}(Bool[
  false,false,true,true,true,false,false,true,
])

### create default events
mod2_events_active_ = NamedTuple{(
  :adm1,:adm2,:adm3,:adm4,
)}(Bool[
  true,true,true,true,
])

### initialization of ODE variables and Records
function mod2_init_func_(cons)
    #(k_a,Vd,Vmax,Km,body_weight,dose_first,dose_last,sigma_add,sigma_prop,ke,) = cons

    # Heta initialize
    t = 0.0 # initial time
    gut = 1.0
    C = 0.0
    BrAC = C / 1.0
    q_a = 0.0
    vabs = cons[1] * q_a
    blood = cons[2] * cons[5]
    vel = cons[10] * C * blood + cons[3] * C / (cons[4] + abs(C)) * blood
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    # save results

    return (
        Float64[
            q_a,
            C * blood,
        ],
        Float64[
            gut,
            blood,
        ]
    )
end

### calculate RHS of ODE
function mod2_ode_func_(du, u, p, t)
    cons = p.constants
    (gut,blood,) = p.static
    (q_a,C_,) = u

    # Heta rules
    C = C_ / blood
    BrAC = C / 1.0
    vabs = cons[1] * q_a
    vel = cons[10] * C * blood + cons[3] * C / (cons[4] + abs(C)) * blood
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    #p.static .= [gut,blood,]
    du .= [
      -vabs,  # dq_a_/dt
      vabs-vel,  # dC_/dt
    ]
end

### output function
function mod2_saving_generator_(outputIds::Vector{Symbol})
    function saving_(u, t, integrator)
        cons = integrator.p.constants
        (gut,blood,) = integrator.p.static
        (q_a,C_,) = u

        # Heta rules
        C = C_ / blood
        BrAC = C / 1.0
        vabs = cons[1] * q_a
        vel = cons[10] * C * blood + cons[3] * C / (cons[4] + abs(C)) * blood
        sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

        # force amount
        q_a_ = q_a

        d = Base.@locals
        return [d[id] for id in outputIds]
    end
end

### TIME EVENTS ###
function mod2_adm1_tstops_func_(cons, times)
    return [0.0]
end
function mod2_adm2_tstops_func_(cons, times)
    return [0.333]
end
function mod2_adm3_tstops_func_(cons, times)
    return [0.667]
end
function mod2_adm4_tstops_func_(cons, times)
    return [1.5]
end

### D EVENTS ###

### STOP EVENTS ###

### event assignments
function mod2_adm1_affect_func_(integrator)
    cons = integrator.p.constants
    t = integrator.t
    (gut,blood,) = integrator.p.static
    (q_a,C_,) = integrator.u

    # Heta rules
    C = C_ / blood
    BrAC = C / 1.0
    vabs = cons[1] * q_a
    vel = cons[10] * C * blood + cons[3] * C / (cons[4] + abs(C)) * blood
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    # recalculated values
    (q_a,) = (q_a + cons[6] * cons[5],)

    # save results
    integrator.u[[1,]] = [
        q_a,
    ]
    integrator.p.static[[]] = [
    ]
end
function mod2_adm2_affect_func_(integrator)
    cons = integrator.p.constants
    t = integrator.t
    (gut,blood,) = integrator.p.static
    (q_a,C_,) = integrator.u

    # Heta rules
    C = C_ / blood
    BrAC = C / 1.0
    vabs = cons[1] * q_a
    vel = cons[10] * C * blood + cons[3] * C / (cons[4] + abs(C)) * blood
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    # recalculated values
    (q_a,) = (q_a + cons[6] * cons[5],)

    # save results
    integrator.u[[1,]] = [
        q_a,
    ]
    integrator.p.static[[]] = [
    ]
end
function mod2_adm3_affect_func_(integrator)
    cons = integrator.p.constants
    t = integrator.t
    (gut,blood,) = integrator.p.static
    (q_a,C_,) = integrator.u

    # Heta rules
    C = C_ / blood
    BrAC = C / 1.0
    vabs = cons[1] * q_a
    vel = cons[10] * C * blood + cons[3] * C / (cons[4] + abs(C)) * blood
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    # recalculated values
    (q_a,) = (q_a + cons[6] * cons[5],)

    # save results
    integrator.u[[1,]] = [
        q_a,
    ]
    integrator.p.static[[]] = [
    ]
end
function mod2_adm4_affect_func_(integrator)
    cons = integrator.p.constants
    t = integrator.t
    (gut,blood,) = integrator.p.static
    (q_a,C_,) = integrator.u

    # Heta rules
    C = C_ / blood
    BrAC = C / 1.0
    vabs = cons[1] * q_a
    vel = cons[10] * C * blood + cons[3] * C / (cons[4] + abs(C)) * blood
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    # recalculated values
    (q_a,) = (q_a + cons[7] * cons[5],)

    # save results
    integrator.u[[1,]] = [
        q_a,
    ]
    integrator.p.static[[]] = [
    ]
end

### MODELS ###

mod2_model_ = (
  mod2_init_func_,
  mod2_ode_func_,
  NamedTuple{(
    :adm1,
    :adm2,
    :adm3,
    :adm4,
  )}([
    (mod2_adm1_tstops_func_, mod2_adm1_affect_func_, :adm1, false),
    (mod2_adm2_tstops_func_, mod2_adm2_affect_func_, :adm2, false),
    (mod2_adm3_tstops_func_, mod2_adm3_affect_func_, :adm3, false),
    (mod2_adm4_tstops_func_, mod2_adm4_affect_func_, :adm4, false),
  ]),
  NamedTuple{(
  )}([
  ]),
  NamedTuple{(
  )}([
  ]),
  mod2_saving_generator_,
  mod2_constants_num_,
  mod2_events_active_,
  mod2_records_output_
)
### create default constants
mod3_constants_num_ = NamedTuple{(
  :k_a,:Vd,:Vmax,:Km,:body_weight,:dose_first,:dose_last,:sigma_add,:sigma_prop,:t_sw,:ke,
)}(Float64[
  0.062,0.457,0.136,0.096,70,0.186,0.112,0.1,0,2.6,0.011,
])

### create default observables
mod3_records_output_ = NamedTuple{(
  :gut,:blood,:q_a,:C,:BrAC,:vabs,:vel,:sigma1,
)}(Bool[
  false,false,true,true,true,false,false,true,
])

### create default events
mod3_events_active_ = NamedTuple{(
  :adm1,:adm2,:adm3,:adm4,
)}(Bool[
  true,true,true,true,
])

### initialization of ODE variables and Records
function mod3_init_func_(cons)
    #(k_a,Vd,Vmax,Km,body_weight,dose_first,dose_last,sigma_add,sigma_prop,t_sw,ke,) = cons

    # Heta initialize
    t = 0.0 # initial time
    gut = 1.0
    C = 0.0
    BrAC = C / 1.0
    q_a = 0.0
    vabs = cons[1] * q_a
    blood = cons[2] * cons[5]
    vel = (t < cons[10]) ? (cons[11] * C * blood) : (cons[3] * C / (cons[4] + abs(C)) * blood)
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    # save results

    return (
        Float64[
            q_a,
            C * blood,
        ],
        Float64[
            gut,
            blood,
        ]
    )
end

### calculate RHS of ODE
function mod3_ode_func_(du, u, p, t)
    cons = p.constants
    (gut,blood,) = p.static
    (q_a,C_,) = u

    # Heta rules
    C = C_ / blood
    BrAC = C / 1.0
    vabs = cons[1] * q_a
    vel = (t < cons[10]) ? (cons[11] * C * blood) : (cons[3] * C / (cons[4] + abs(C)) * blood)
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    #p.static .= [gut,blood,]
    du .= [
      -vabs,  # dq_a_/dt
      vabs-vel,  # dC_/dt
    ]
end

### output function
function mod3_saving_generator_(outputIds::Vector{Symbol})
    function saving_(u, t, integrator)
        cons = integrator.p.constants
        (gut,blood,) = integrator.p.static
        (q_a,C_,) = u

        # Heta rules
        C = C_ / blood
        BrAC = C / 1.0
        vabs = cons[1] * q_a
        vel = (t < cons[10]) ? (cons[11] * C * blood) : (cons[3] * C / (cons[4] + abs(C)) * blood)
        sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

        # force amount
        q_a_ = q_a

        d = Base.@locals
        return [d[id] for id in outputIds]
    end
end

### TIME EVENTS ###
function mod3_adm1_tstops_func_(cons, times)
    return [0.0]
end
function mod3_adm2_tstops_func_(cons, times)
    return [0.333]
end
function mod3_adm3_tstops_func_(cons, times)
    return [0.667]
end
function mod3_adm4_tstops_func_(cons, times)
    return [1.5]
end

### D EVENTS ###

### STOP EVENTS ###

### event assignments
function mod3_adm1_affect_func_(integrator)
    cons = integrator.p.constants
    t = integrator.t
    (gut,blood,) = integrator.p.static
    (q_a,C_,) = integrator.u

    # Heta rules
    C = C_ / blood
    BrAC = C / 1.0
    vabs = cons[1] * q_a
    vel = (t < cons[10]) ? (cons[11] * C * blood) : (cons[3] * C / (cons[4] + abs(C)) * blood)
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    # recalculated values
    (q_a,) = (q_a + cons[6] * cons[5],)

    # save results
    integrator.u[[1,]] = [
        q_a,
    ]
    integrator.p.static[[]] = [
    ]
end
function mod3_adm2_affect_func_(integrator)
    cons = integrator.p.constants
    t = integrator.t
    (gut,blood,) = integrator.p.static
    (q_a,C_,) = integrator.u

    # Heta rules
    C = C_ / blood
    BrAC = C / 1.0
    vabs = cons[1] * q_a
    vel = (t < cons[10]) ? (cons[11] * C * blood) : (cons[3] * C / (cons[4] + abs(C)) * blood)
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    # recalculated values
    (q_a,) = (q_a + cons[6] * cons[5],)

    # save results
    integrator.u[[1,]] = [
        q_a,
    ]
    integrator.p.static[[]] = [
    ]
end
function mod3_adm3_affect_func_(integrator)
    cons = integrator.p.constants
    t = integrator.t
    (gut,blood,) = integrator.p.static
    (q_a,C_,) = integrator.u

    # Heta rules
    C = C_ / blood
    BrAC = C / 1.0
    vabs = cons[1] * q_a
    vel = (t < cons[10]) ? (cons[11] * C * blood) : (cons[3] * C / (cons[4] + abs(C)) * blood)
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    # recalculated values
    (q_a,) = (q_a + cons[6] * cons[5],)

    # save results
    integrator.u[[1,]] = [
        q_a,
    ]
    integrator.p.static[[]] = [
    ]
end
function mod3_adm4_affect_func_(integrator)
    cons = integrator.p.constants
    t = integrator.t
    (gut,blood,) = integrator.p.static
    (q_a,C_,) = integrator.u

    # Heta rules
    C = C_ / blood
    BrAC = C / 1.0
    vabs = cons[1] * q_a
    vel = (t < cons[10]) ? (cons[11] * C * blood) : (cons[3] * C / (cons[4] + abs(C)) * blood)
    sigma1 = sqrt(cons[8] ^ 2.0 + C ^ 2.0 * cons[9] ^ 2.0)

    # recalculated values
    (q_a,) = (q_a + cons[7] * cons[5],)

    # save results
    integrator.u[[1,]] = [
        q_a,
    ]
    integrator.p.static[[]] = [
    ]
end

### MODELS ###

mod3_model_ = (
  mod3_init_func_,
  mod3_ode_func_,
  NamedTuple{(
    :adm1,
    :adm2,
    :adm3,
    :adm4,
  )}([
    (mod3_adm1_tstops_func_, mod3_adm1_affect_func_, :adm1, false),
    (mod3_adm2_tstops_func_, mod3_adm2_affect_func_, :adm2, false),
    (mod3_adm3_tstops_func_, mod3_adm3_affect_func_, :adm3, false),
    (mod3_adm4_tstops_func_, mod3_adm4_affect_func_, :adm4, false),
  ]),
  NamedTuple{(
  )}([
  ]),
  NamedTuple{(
  )}([
  ]),
  mod3_saving_generator_,
  mod3_constants_num_,
  mod3_events_active_,
  mod3_records_output_
)

### OUTPUT ###

return (
  (
    nameless = nameless_model_,
    mod1 = mod1_model_,
    mod2 = mod2_model_,
    mod3 = mod3_model_,
  ),
  (),
  "0.6.6"
)

end
)()

max_stars_repo_name: ivborissov/alco

max_issues_repo_name: ivborissov/alco

max_forks_repo_name: ivborissov/alco


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "cilkrts"
version = v"0.1.2"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/OpenCilk/cilkrts.git", "507a93234a480431ee8f7e49bbb464e3572618f5")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/cilkrts/
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN%.*}_gcc.cmake -DCMAKE_BUILD_TYPE=Release ..
make
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("i686", "linux"; libc = "glibc"),
    Platform("x86_64", "linux"; libc = "glibc"),
    Platform("aarch64", "linux"; libc = "glibc"),
    Platform("armv7l", "linux"; call_abi = "eabihf", libc = "glibc"),
    Platform("powerpc64le", "linux"; libc = "glibc"),
    Platform("x86_64", "macos"),
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libcilkrts", :libcilkrts)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"7.1.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@testset "$TEST $G" begin

if !@isdefined(test_find_path_types)
    # Test the types of the values returned by fetch_path
    function test_find_path_types(residual_graph, s, t, flow_matrix, capacity_matrix)
        v, P, S, flag = Erdos.fetch_path(residual_graph, s, t, flow_matrix, capacity_matrix)
        @test typeof(P) == Vector{Int}
        @test typeof(S) == Vector{Int}
        @test typeof(flag) == Int
        @test typeof(v) <: Integer
    end

    # Test the value of the flags returned.
    function test_find_path_disconnected(residual_graph, s, t, flow_matrix, capacity_matrix)
        h = copy(residual_graph)
        for dst in collect(neighbors(residual_graph, s))
            rem_edge!(residual_graph, s, dst)
        end
        v, P, S, flag = Erdos.fetch_path(residual_graph, s, t, flow_matrix, capacity_matrix)
        @test flag == 1
        for dst in collect(neighbors(h, t))
            rem_edge!(h, t, dst)
        end
        v, P, S, flag = Erdos.fetch_path(h, s, t, flow_matrix, capacity_matrix)
        @test flag == 0
        for i in collect(in_neighbors(h, t))
            rem_edge!(h, i, t)
        end
        v, P, S, flag = Erdos.fetch_path(h, s, t, flow_matrix, capacity_matrix)
        @test flag == 2
    end
end

# Construct DiGraph
flow_graph = DG(8)

# Load custom dataset
flow_edges = [
    (1,2,10),(1,3,5),(1,4,15),(2,3,4),(2,5,9),
    (2,6,15),(3,4,4),(3,6,8),(4,7,16),(5,6,15),
    (5,8,10),(6,7,15),(6,8,10),(7,3,6),(7,8,10)
]

capacity_matrix = zeros(Int,8,8)

for e in flow_edges
    u,v,f = e
    add_edge!(flow_graph,u,v)
    capacity_matrix[u,v] = f
end

residual_graph = complete(flow_graph)

# Test with default distances
@test Erdos.edmonds_karp_impl(residual_graph, 1, 8, Erdos.DefaultCapacity(residual_graph))[1] == 3

# Test with capacity matrix
@test Erdos.edmonds_karp_impl(residual_graph,1,8,capacity_matrix)[1] == 28

flow_matrix = zeros(Int, nv(residual_graph), nv(residual_graph))
test_find_path_types(residual_graph, 1,8, flow_matrix, capacity_matrix)
test_find_path_disconnected(residual_graph, 1, 8, flow_matrix, capacity_matrix)

end # testset


text:
formatter_options = (
  indent = 2,
  margin = 80,
  always_for_in = true,
  whitespace_typedefs = true,
  whitespace_ops_in_indices = true,
  remove_extra_newlines = true,
)


text:
using PFFRGSolver

println()

# test frequency implementation
println("Running frequency tests ...")
test_frequencies()
println()

# test action implementation
println("Running lattice tests ...")
test_lattice()
println()

# test action implementation
println("Running action tests ...")
test_action()
println()

# test flow implementation
println("Running flow tests ...")
test_flow()
println()

# test observable implementation
println("Running observable tests ...")
test_observable()
println()


text:
using Metaheuristics
using Test
import Random: seed!

@testset "Common Methods" begin
    function simple_test()

        # compare
        x = rand(10)

        @test Metaheuristics.compare(x, x) == 0
        @test Metaheuristics.compare(x, x .+ 1) == 1
        @test Metaheuristics.compare(x .+ 1, x) == 2
        y = copy(x)
        y[1] -= 1.0
        x[end] -= 1.0
        @test Metaheuristics.compare(x, y) == 3

        sol_feasible   = Metaheuristics.create_child(x, (x, zeros(2), zeros(2)))
        sol_infeasible = Metaheuristics.create_child(x, (x, rand(4), rand(8)))

        @test Metaheuristics.compare(sol_feasible, sol_infeasible) == 1
        @test Metaheuristics.compare(sol_infeasible, sol_feasible) == 2
        @test Metaheuristics.is_better(sol_feasible, sol_infeasible) # true
        @test !Metaheuristics.is_better(sol_infeasible, sol_feasible)# !false == true

        # repair
        x = collect(range(-10, 10, length=20))
        bounds = [-ones(20) ones(20)]'
        for repair in [Metaheuristics.reflected_back_to_bound!,
                       Metaheuristics.replace_with_random_in_bounds!,
                       Metaheuristics.wrap_to_bounds!,
                       Metaheuristics.reset_to_violated_bounds!,
                      ]
            @test Metaheuristics.is_in_bounds(repair(x, bounds), bounds)
        end
        @test Metaheuristics.is_in_bounds(Metaheuristics.evo_boundary_repairer!(x, zeros(20), bounds), bounds)

    end

    function test_problems()

        # single-objective
        for problem = [:sphere, :discus, :rastrigin, ]
            f, bounds, optimums = Metaheuristics.TestProblems.get_problem(problem)
            for opt in optimums
                @test f(opt.x) == opt.f == 0.0
            end
        end

        # single-objective constrained
        for problem = [:constrained1, :constrained2, :constrained3]
            f, bounds, optimums = Metaheuristics.TestProblems.get_problem(problem)
            for opt in optimums
                fx, gx, hx = f(opt.x)
                @test fx == opt.f && sum(max.(gx, 0.0)) ≈ 0.0 && sum(abs.(hx)) ≈ 0.0
            end
        end

        # multi-objective
        for problem = [:ZDT1, :ZDT2, :ZDT3, :ZDT4, :ZDT6,
                       :DTLZ1, :DTLZ2, :DTLZ3, :DTLZ4, :DTLZ5, :DTLZ6]
            f, bounds, optimums = Metaheuristics.TestProblems.get_problem(problem)
            pf = Metaheuristics.get_non_dominated_solutions(optimums)

            @test length(pf) == length(optimums)

            # test dimensions
            fx, gx, hx = f(bounds[1,:])
            @test length(fx) == length(optimums[1].f)

            # test performance indicators
            ## generational distance
            @test Metaheuristics.PerformanceIndicators.igd(optimums, pf) ≈ 0.0
            @test Metaheuristics.PerformanceIndicators.igd_plus(optimums, pf) ≈ 0.0
            @test Metaheuristics.PerformanceIndicators.gd(optimums, pf) ≈ 0.0
            @test Metaheuristics.PerformanceIndicators.gd_plus(optimums, pf) ≈ 0.0

            ## spacing
            seed!(1)
            X = rand(100, 7)

            @test Metaheuristics.PerformanceIndicators.spacing(X) < 0.5
            @test Metaheuristics.PerformanceIndicators.spacing(optimums) < 0.2

            ## covering
            ## non dominated solutions
            X = fvals(pf)
            Y = copy(X)
            # 50% of solutions in Y will be dominated by X
            mask = 1:size(X,1) ÷ 2
            Y[mask,1] = Y[mask,1] .+ 1.0
            Cxy=Metaheuristics.PerformanceIndicators.covering(X, Y)
            Cyx=Metaheuristics.PerformanceIndicators.covering(Y, X)
            @test Cxy == (size(X,1) ÷ 2) / size(X,1)
            @test Cyx == 0.0

            ## Δₚ (delta p)
            @test PerformanceIndicators.deltap(pf, optimums) == PerformanceIndicators.Δₚ(pf, optimums)
        end

    end

    function test_hypervolume()
        front = [[0.        ,0.        ,1.1       ],
                 [0.        ,1.07863874,0.21572775],
                 [0.13540064,0.94780448,0.54160256],
                 [0.29938208,0.7484552 ,0.7484552 ],
                 [0.46669048,0.62225397,0.77781746],
                 [0.6350853 ,0.6350853 ,0.6350853 ],
                 [0.7484552 ,0.7484552 ,0.29938208],
                 [0.89510682,0.        ,0.63936201],
                 [0.98386991,0.49193496,0.        ],
                 [2.0, 3, 2], # only for testing
                 [1.1       ,0.        ,0.        ]]

        referencePoint = [2.0, 2, 2]
        hyperVolume = Metaheuristics.PerformanceIndicators.hypervolume(front, referencePoint)
        @test hyperVolume ≈ 6.793879034744429

        front_ = Array(hcat(front...)')
        hyperVolume = Metaheuristics.PerformanceIndicators.hypervolume(front_, referencePoint)
        @test hyperVolume ≈ 6.793879034744429

        # testing nadir and ideal pint
        @test ideal(front) == ideal(Array(hcat(front...)')) == zeros(3)
        @test nadir(front) == nadir(Array(hcat(front...)')) == [2.0, 3, 2]

        f, bounds, front2 = Metaheuristics.TestProblems.get_problem(:DTLZ2);
        @test Metaheuristics.PerformanceIndicators.hypervolume(front2, nadir(front2)) > 0
        @test ideal(front2) != nadir(front2)
    end

    function test_prev_pop()
        f, bounds, optimums = Metaheuristics.TestProblems.get_problem(:sphere)

        options=Options(f_calls_limit=10000, seed=1)

        method = DE(options = options)
        res = optimize(f, bounds, method)
        optimum = minimum(res)
        N = length(res.population)

        methods = [ECA(N=N,options = options),
                   PSO(N=N,options = options),
                   DE(N=N,options = options)]

        for method in methods
            prev_status = State(res.best_sol, res.population)
            method.status = prev_status
            res = optimize(f, bounds, method)
            optimum = minimum(res)
        end

        @test optimum < 1e-5
    end

    simple_test()
    test_problems()
    test_prev_pop()
    test_hypervolume()

end


text:
"""
# Usage

    solve(prob::MultiTermsFODEProblem, h, FODEMatrixDiscrete())

Using [triangular strip matrices](https://en.wikipedia.org/wiki/Triangular_matrix) to discrete fractional ordinary differential equations to simple algebra system and solve the system.

### References

```tex
@inproceedings{Podlubny2000MATRIXAT,
  title={MATRIX APPROACH TO DISCRETE FRACTIONAL CALCULUS},
  author={Igor Podlubny},
  year={2000}
}
```
"""
struct FODEMatrixDiscrete <: FractionalDiffEqAlgorithm end

"""
# Usage

    solve(α, β, κ, T, M, N, FPDEMatrixDiscrete())

Using [triangular strip matrices](https://en.wikipedia.org/wiki/Triangular_matrix) to discrete fractional partial differential equations to simple algebra system and solve the system.

```tex
@article{2009,
   title={Matrix approach to discrete fractional calculus II: Partial fractional differential equations},
   DOI={10.1016/j.jcp.2009.01.014},
   author={Podlubny, Igor and Chechkin, Aleksei and Skovranek, Tomas and Chen, YangQuan and Vinagre Jara, Blas M.},
}
```
"""
struct FPDEMatrixDiscrete <: FractionalDiffEqAlgorithm end

function solve(prob::MultiTermsFODEProblem, h, ::FODEMatrixDiscrete)
    @unpack parameters, orders, rightfun, tspan = prob
    t0 = tspan[1]; T = tspan[2]
    N::Int64 = floor(Int, (T-t0)/h)
    highestorder = Int64(findmax(ceil.(orders))[1])
    rows = collect(Int64, 1:highestorder)

    equation = zeros(N, N)

    for (i, j) in zip(parameters, orders)
        equation += i*D(N, j, h)
    end

    equation = eliminator(N, rows)*equation*eliminator(N, rows)'

    if typeof(rightfun) <: Number # Handling right hand side
        rightside = eliminator(N, rows)*rightfun*ones(N)
    else
        rightside = eliminator(N, rows)*(rightfun.(collect(h:h:T)).+ic_handling(orders, parameters, u0))
    end

    result = equation\rightside
    result = vcat(zeros(highestorder), result)

    return FODESolution(collect(Float64, t0+h:h:T), result)
end

"""
    eliminator(n, row)

Compute the eliminator matrix Sₖ by omiting n-th row
"""
function eliminator(n, row)
    temp = zeros(n, n) + I
    return temp[Not(row), :]
end

"""
Generating elements in Triangular Strip Matrix.
"""
function omega(n, α)
    omega = zeros(n+1)

    omega[1] = 1
    @fastmath @inbounds @simd for i ∈ 1:n
        omega[i+1] = (1 - (α+1)/i)*omega[i]
    end

    return omega
end

"""
    D(N, α, h)

Using D function to construct left hand side equations.

!!! info
    Here ```N``` is the size of discrete matrix.
"""
function D(N, α, h)
    α == 0 ? (return zeros(N, N) + I) : nothing # When α=0, D is an identity matrix.
    result = zeros(N, N)
    temp = omega(N, α)

    @fastmath @inbounds @simd for i in range(1, N, step=1)
        result[i, 1:i] = reverse(temp[1:i])
    end

    return h^(-α)*result
end

function F(N, α, h)
    result = zeros(N, N)
    temp = omega(N, α)

    @fastmath @inbounds @simd for i ∈ 1:N
        result[i, 1:i] = reverse(temp[1:i])
    end

    result = reverse(reverse(result, dims=1), dims=2)

    return (-1)^(ceil(α))*h^(-α)*result
end

function meshgrid(xin,yin)
    nx = length(xin)
    ny = length(yin)
    xout = zeros(ny,nx)
    yout = zeros(ny,nx)
    for jx = 1:nx
        for ix = 1:ny
            xout[ix,jx] = xin[jx]
            yout[ix,jx] = yin[ix]
        end
    end
    return (x=xout, y=yout)
end

# Construct Riesz Symmetric Matrix
function RieszMatrix(α, N, h)
    caputo = B(N+1, α)
    caputo = caputo[2:(N+1), 1:N]
    result = 1/2*(caputo+caputo')
    result = h^(-α)*result

    return result
end

function B(N, p)
    result = zeros(N, N)
    temp = omega(N, p)

    @inbounds @simd for i ∈ 1:N
        @views result[i, 1:i]=reverse(temp[1:i])
    end

    return result
end

"""
    bagleytorvik(p1, p2, p3, T, h)

By specifying the parameters of Bagley Torvik Equation, we can use **bagleytorvik** to directly obtain the numerical approximation of a bagley torvik equation.

!!! info "p2 ≠ 0"
    Please note that the parameter of fractional derivative part must not be 0
"""
function bagleytorvik(p1, p2, p3, right, T, h)
    N = round(Int, T/h)+1
    equation = p1*D(N, 2, h) + p2*D(N, 1.5, h) + p3*(zeros(N, N)+I)
    equation = eliminator(N, [1,2])*equation*eliminator(N, [1,2])'

    rows = collect(Int64, 1:2)

    typeof(right) <: Number ? rightside = eliminator(N, rows)*right*ones(N) : rightside = eliminator(N, rows)*right.(collect(h:h:T))

    result = equation\rightside

    return vcat(zeros(Float64, 2), result)
end

function solve(α, β, κ, T, M, N, ::FPDEMatrixDiscrete)
    h = T/(M-1)
    τ = h^2/6

    # Construct the time partial derivative matrix
    TMatrix = kron(D(N-1, α, τ)', zeros(M, M) + I)

    # Construct the spatial partial derivative matrix
    SMatrix = kron(zeros(N-1, N-1) + I, RieszMatrix(β, M, h))

    system = TMatrix - κ.*SMatrix

    # Handling boundary conditions
    BMatrix = kron(zeros(N-1, N-1) + I, eliminator(M, [1, M]))
    system = system*BMatrix'

    left = BMatrix*system

    tmp = left\(ones(size(left, 1), 1))

    YS = reshape(tmp, M-2, N-1)
    YS = reverse(YS, dims=2)
    U = YS
    rows, columns = size(U)

    # Boundry conditions
    U = [zeros(1, columns); U; zeros(1, columns)]
    result = [zeros(M) U]

    return result
end
#=
"""
Shifter matrices for representing time delays with discrete matrix.
"""
function shift(U, k)
    s = size(U)
    n = s[1]
    absk = abs(k)
    ek = diagm(k => ones(n-absk))
    s1 = eliminator(n, collect(1:absk))
    s2 = eliminator(n, collect(n-absk+1:n))

    if k>0
        Y = s2*ek*U*ek*s1'
    elseif k<0
        Y = s1*ek*U*ek*s2'
    else
        Y = U
    end
    return Y
end

#FIXME: Fractional partial differential equations with time delay.
function fdpde(alpha, alphad, beta, steps)
    a2=1
    L = 1
    m = 21
    n = 148
    h = L/(m-1)
    tau = h^2 / (6*a2)

    B1 = D(n-1, alpha, tau)'
    TD = kron(B1, zeros(m, m) + I)

    Bdelay = shift(D(n-1+steps, alphad, tau)', steps)
    TDdelay = kron(Bdelay, zeros(m, m) + I)

    B2 = RieszMatrix(beta, m, h)
    SD = kron(zeros(n-1, n-1) + I, B2)

    SystemMatrix = 0.5*TD + 0.5*TDdelay- a2*SD

    S = eliminator(m, [1, m])
    SK = kron(zeros(n-1, n-1) + I, S)
    SystemMatrix_without_columns_1_m = SystemMatrix * SK'

    S = eliminator(m, [1 m])
    SK = kron(zeros(n-1, n-1) + I, S)
    SystemMatrix_without_rows_columns_1_m = SK * SystemMatrix_without_columns_1_m

    F = 8*ones(size(SystemMatrix_without_rows_columns_1_m,1), 1)

    Y = SystemMatrix_without_rows_columns_1_m\F

    YS = reshape(Y, m-2, n-1)
    YS = reverse(YS, dims=2)

    U = YS

    (rows, columns) = size(U)
    U = [zeros(1, columns); U;  zeros(1, columns)]
    U = [zeros(1,m)' U];
    (XX, YY)=meshgrid(tau*(0:n-1),h*(0:m-1))
    return XX, YY, U
end
=#
#=
sol=fdpde(0.5, 0.5, 0.5, 1)
using Plots
plotlyjs()
plot(sol, st=:surface)
=#

#=
Test code for FPDEMatrixDiscrete, to explore the wide usage of FPDEMatrixDiscrete
=#
#=
function testsolve(α, β, T, M, N)
    h = T/(M-1)
    τ = h^2/6

    # Construct the time partial derivative matrix
    TMatrix = kron(D(N-1, α, τ)', zeros(M, M) + I)

    # Construct the spatial partial derivative matrix
    SMatrix = kron(zeros(N-1, N-1) + I, RieszMatrix(β, M, h))

    system = TMatrix+SMatrix

    #FIXME: Handling boundary conditions
    BMatrix = kron(zeros(N-1, N-1) + I, eliminator(M, [1, M]))
    system = system*BMatrix'

    leftside = BMatrix*system

    XX = range(1/2793, 1, step=1/2793)
    YY = range(1/2793, 1, step=1/2793)

    result = leftside\(-sin.(π .*XX).*sin.(π .*YY))

    return result
end

using Plots

using LinearAlgebra, InvertedIndices

tmp = testsolve(2, 2, 1, 21, 148)

YS = reshape(tmp, 19, 147)
YS = reverse(YS, dims=2)
U = YS

rows, columns = size(U)

U = [zeros(1, columns); U; zeros(1, columns)]
U=[zeros(1, 21)' U]

XX, YY = meshgrid(0.05^2/6 .*(0:147), 0:0.05:1)

plotlyjs()

plot(XX, YY, U, st=:surface)

=#

"""
    diffusion(α, β)

**Diffusion equation** with time fractional derivative.

α and β are the coefficients of the time derivative and spatial derivative.
"""
function diffusion(α, β)
    solve(α, β, 1, 21, 148, FPDEMatrixDiscrete())
end

## An fractional partial differential equation Example


text:
# Sample postprocessing script for virtual satellite tracking.
# Outputs are stored in CSV format for sharing data.
#
# Usage:
#   julia -t nthreads demo_virtual_satellite.jl
# or
#   JULIA_NUM_THREADS=nthreads julia demo_virtual_satellite.jl
#
# Hongyang Zhou, hyzhou@umich.edu

using Glob, DelimitedFiles, Vlasiator, DataFrames

function extract_vars(files, loc)
   nfiles = length(files)
   # variables to be extracted
   t   = zeros(Float32, nfiles)
   rho = zeros(Float32, nfiles)
   vx  = zeros(Float32, nfiles)
   vy  = zeros(Float32, nfiles)
   p   = zeros(Float32, nfiles)
   bz  = zeros(Float32, nfiles)
   ex  = zeros(Float32, nfiles)
   ey  = zeros(Float32, nfiles)

   local id
   let meta = load(files[1])
      id = getcell(meta, loc)
   end

   # Extract data from each frame
   Threads.@threads for i = eachindex(files)
      meta = load(files[i])
      t[i] = meta.time
      rho[i] = readvariable(meta, "proton/vg_rho", id)[1]
      vx[i], vy[i] = readvariable(meta, "proton/vg_v", id)[1:2]
      p[i] = readvariable(meta, "vg_pressure", id)[1]
      bz[i] = readvariable(meta, "vg_b_vol", id)[3]
      ex[i], ey[i] = readvariable(meta, "vg_e_vol", id)[1:2]
   end

   df = DataFrame(t = t, rho = rho, vx = vx, vy = vy, p = p, bz = bz, ex = ex, ey = ey)
   # Save into text file
   writedlm("satellite.csv", Iterators.flatten(([names(df)], eachrow(df))), ',')
end

#####

Re = Vlasiator.Re # Earth radius

files = glob("bulk*.vlsv", "./")

# virtual satellite location
loc = [12Re, 0, 0]

println("Number of files: $(length(files))")
println("Extracting location: $loc")
println("Running with $(Threads.nthreads()) threads...")

@time extract_vars(files, loc)

println("Virtual satellite extraction done!")

## Visualization
#=
using PyPlot, DelimitedFiles

data = readdlm("satellite.csv", ','; header=true)

fig, ax = subplots(figsize=(8,10), 5,1, sharex=true, constrained_layout=true)

ax[1].plot(data[1][:,1], data[1][:,2] ./ 1e6, label="density")
ax[2].plot(data[1][:,1], data[1][:,3] ./ 1e3, label="vx")
ax[2].plot(data[1][:,1], data[1][:,4] ./ 1e3, label="vy")
ax[3].plot(data[1][:,1], data[1][:,5] .* 1e9, label="p")
ax[4].plot(data[1][:,1], data[1][:,6] .* 1e9, label="bz")
ax[5].plot(data[1][:,1], data[1][:,7] .* 1e3, label="ex")
ax[5].plot(data[1][:,1], data[1][:,8] .* 1e3, label="ey")

for a in ax
   a.grid(true)
   a.legend()
end

ax[1].set_ylabel("density [amu/cc]", fontsize=14)
ax[2].set_ylabel("velocity [km/s]", fontsize=14)
ax[3].set_ylabel("pressure [nPa]", fontsize=14)
ax[4].set_ylabel("magnetic field [nT]", fontsize=14)
ax[5].set_ylabel("electric field [mV/m]", fontsize=14)
ax[5].set_xlabel("time [s]", fontsize=14)

fig.suptitle("Density Pulse Run, location = [12, 0, 0]", fontsize="xx-large")

savefig("virtual_satellite.png", bbox_inches="tight")
=#