text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

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

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

@inline function exp_ps_m6_opt_rho2_22!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=9
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 T2k9
    # Computing Bb2 = x*I+x*A
    coeff1=0.6932267294704894
    coeff2=0.42019701808349896
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.6932267294704894
    coeff2=0.42019701808349896
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=0.14883975012114015
    coeff2=0.8408521785457064
    coeff3=0.6441812976652185
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=1.042613839944681
    coeff2=0.3549827751063252
    coeff3=0.19349580056392088
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.03430490914845939
    coeff2=0.3596800085189574
    coeff3=0.8532569213076548
    coeff4=0.9021858527038544
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=1.0264634378476645
    coeff2=0.6572335936947115
    coeff3=-0.054585888069760787
    coeff4=0.028297918922423518
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=1.0499238030144577e-10
    coeff2=0.0012477776392267876
    coeff3=0.000784769897476988
    coeff4=-0.10903345552815719
    coeff5=0.9996647988764669
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=4.054924703391629e-7
    coeff2=-2.3045432301617836e-8
    coeff3=-7.77392066692997e-10
    coeff4=3.2001904605874257e-10
    coeff5=-1.1346120403005061e-12
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=2.880470058346373e-5
    coeff2=-0.04253750167360615
    coeff3=-0.13303235916009215
    coeff4=0.8152592432110274
    coeff5=0.9444261197925446
    coeff6=1.7522893501377047e-7
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.0004535421655560562
    coeff2=2.887116538841799e-5
    coeff3=2.7169875443715466e-5
    coeff4=-1.1358540783808526e-6
    coeff5=1.232855953231344e-7
    coeff6=1.0000000017458077
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing Ba7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.9986923999012539
    coeff2=0.04556161788104884
    coeff3=0.39768438765561237
    coeff4=0.22180315416394053
    coeff5=1.094529320861645
    coeff6=-1.9466631013553095e-7
    coeff7=0.00024920094919997024
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=-0.03834000555185737
    coeff2=0.014143811582239133
    coeff3=0.008028642819640242
    coeff4=0.0023561023048899606
    coeff5=3.575682812162475e-5
    coeff6=4.24521104018061e-5
    coeff7=1.0000015155621724
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots9,memslots2,memslots3)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 3
    # Computing T2k9 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.7426007763589999
    coeff2=0.23958838183794046
    coeff3=0.2456652487277263
    coeff4=0.19442822324865222
    coeff5=0.12981790073255753
    coeff6=2.3660170281050782e-7
    coeff7=0.0013719348524986168
    coeff8=0.9988773313640329
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    # Deallocating B7 in slot 9
    return memslots1 # Returning T2k9
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

fmi3Load("BouncingBall", "ModelicaReferenceFMUs", "0.0.14")

inst = fmi3InstantiateCoSimulation!(myFMU; loggingOn=false)
@test inst != 0

@test fmi3EnterInitializationMode(inst) == 0
@test fmi3ExitInitializationMode(inst) == 0

numStates = length(myFMU.modelDescription.stateValueReferences)
targetValues = [[0.0, 0.0], [1.0, 0.0]]
dir_ders_buffer = zeros(fmi3Float64, numStates)
sample_ders_buffer = zeros(fmi3Float64, numStates, 1)
for i in 1:numStates

    if fmi3ProvidesDirectionalDerivatives(myFMU)
        # multi derivatives calls
        sample_ders = fmi3SampleDirectionalDerivative(inst, myFMU.modelDescription.derivativeValueReferences, [myFMU.modelDescription.stateValueReferences[i]])
        fmi3SampleDirectionalDerivative!(inst, myFMU.modelDescription.derivativeValueReferences, [myFMU.modelDescription.stateValueReferences[i]], sample_ders_buffer)

        @test sum(abs.(sample_ders[:,1] - targetValues[i])) < 1e-3
        @test sum(abs.(sample_ders_buffer[:,1] - targetValues[i])) < 1e-3

        dir_ders = fmi3GetDirectionalDerivative(inst, myFMU.modelDescription.derivativeValueReferences, [myFMU.modelDescription.stateValueReferences[i]])
        fmi3GetDirectionalDerivative!(inst, myFMU.modelDescription.derivativeValueReferences, [myFMU.modelDescription.stateValueReferences[i]], dir_ders_buffer)

        @test sum(abs.(dir_ders - targetValues[i])) < 1e-3
        @test sum(abs.(dir_ders_buffer - targetValues[i])) < 1e-3

        # single derivative call
        dir_der = fmi3GetDirectionalDerivative(inst, myFMU.modelDescription.derivativeValueReferences[1], myFMU.modelDescription.stateValueReferences[1])
        @test dir_der == targetValues[1][1]
    else
        @warn "Skipping directional derivative testing, this FMU from $(ENV["EXPORTINGTOOL"]) doesn't support directional derivatives."
    end
end

# Bug in the FMU
# jac = fmi3GetJacobian(inst, myFMU.modelDescription.derivativeValueReferences, myFMU.modelDescription.stateValueReferences)
# @test jac ≈ hcat(targetValues...)

# jac = fmi3SampleDirectionalDerivative(inst, myFMU.modelDescription.derivativeValueReferences, myFMU.modelDescription.stateValueReferences)
# @test jac ≈ hcat(targetValues...)

fmi3Unload(myFMU)

ext: jl


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -2.9521310086708863 <= q <= 4.804083818674538)

                     add_NL_constraint(m, :(softplus(0.9022086576378046 + 0.4902152360110663*$(x[1]) + 0.6884954091678397*$(x[2]) + 0.4475675674383903*$(x[3])) + softplus(0.023767747364021208 + -0.6986395659548528*$(x[1]) + 0.5778907575747945*$(x[2]) + 0.9752988775257685*$(x[3])) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function formatinfo(s::Spec, x::BigFloat)
    str = Base.MPFR.string_mpfr(abs(s.type == '%' ? 100x : x), makefmt(s))
    n = ncodeunits(str)
    m = 0  # the number of integral digits
    while m + 1 ≤ n && codeunit(str, m + 1) != UInt8('.')
        m += 1
    end
    width = n
    if s.sign == SIGN_PLUS || s.sign == SIGN_SPACE || s.sign != SIGN_NONE && signbit(x)
        width += 1
    end
    if isspecified(s.grouping)
        width += div(m - 1, 3)  # separators
    end
    if !isspecified(s.type) && isinteger(x)
        width += 2
    end
    return width + paddingsize(s, width), str
end

function formatfield(data::Vector{UInt8}, p::Int, s::Spec, x::BigFloat, str::String)
    start = p
    p, signed = sign(data, p, x, s.sign)
    n = ncodeunits(str)
    copyto!(data, p, codeunits(str), 1, n)
    p += n
    if !isspecified(s.type) && all(isdigit, str)
        data[p] = UInt8('.')
        p += 1
    end
    return processfloat(data, start, p, s, signed)
end

function makefmt(s::Spec)
    if s.type == 'G' || s.type == 'g' || !isspecified(s.type)
        type = default(s.type, 'g')
        return isspecified(s.precision) ? "%.$(s.precision)R$(type)" : "%R$(type)"
    elseif s.type == '%'
        precision = default(s.precision, 6)
        return "%.$(precision)Rf"
    else
        precision = default(s.precision, 6)
        return "%.$(precision)R$(s.type)"
    end
end


text:
@testset "NumField/Elem" begin
  Qx, x = PolynomialRing(FlintQQ, "x")
  QasNumberField, _ = NumberField(x - 1)
  Kt, t = PolynomialRing(QasNumberField, "t")
  K1, a1 = NumberField(x^3 - 2)
  K2, (a2, ) = NumberField([x^3 - 2])
  K3, a3 = NumberField(t^3 - 2)
  K4, (a4, ) = NumberField([t^3 - 2])

  @testset for (K, a) in [(K1, a1), (K2, a2), (K3, a3), (K4, a4)]
    b = one(K)
    fl = @inferred isintegral(b)
    @test fl
    fl = @inferred isintegral(a)
    OK = maximal_order(K)
    @test isintegral(a)
    @test isintegral(b)

    B = @inferred basis(K)
    c = @inferred basis_matrix([one(K), a^4])
    @assert nrows(c) == 2
    @assert ncols(c) == 3
    @assert sum(B[i] * c[1, i] for i in 1:3) == one(K)
    @assert sum(B[i] * c[2, i] for i in 1:3) == a^4

    f = @inferred charpoly(one(K()))
    @test f == (gen(parent(f)) - 1)^3
    f = @inferred minpoly(one(K))
    @test f == (gen(parent(f)) - 1)

    f = @inferred charpoly(a)
    @test f == gen(parent(f))^3 - 2
    f = @inferred minpoly(a)
    @test f == gen(parent(f))^3 - 2

  end

  @testset "NumField/Component" begin
    Qx, x = PolynomialRing(FlintQQ, "x")
    QasNumberField, _ = NumberField(x - 1)
    Kt, t = PolynomialRing(QasNumberField, "t")
    K3, a3 = NumberField(t^3 - 2)
    K4, (a4, ) = NumberField([t^3 - 2])
    K41, mK4 = component(K4, 1)
    @test isisomorphic(K41, K3)[1]
    @test mK4(gen(K41)) == a4

    K1, a1 = NumberField(x^3 - 2)
    K2, (a2, ) = NumberField([x^3 - 2])
    K21, mK2 = component(K2, 1)
    @test isisomorphic(K21, K1)[1]
    @test mK2(gen(K21)) == a2
  end

  @testset "rand" begin
    Qx, x = PolynomialRing(FlintQQ, "x")
    K, a = NumberField(x^3 - 2)
    v = [a^0, a^2]
    @assert elem_type(K) == nf_elem
    for args = ((v, 1:3), (v, 1:3, 2))
      m = make(K, args...)
      for x in (rand(args...), rand(rng, args...),
                rand(m), rand(rng, m))
        @test x isa nf_elem
      end
      @test rand(m, 3) isa Vector{nf_elem}
      c = zero(K)
      @test c === rand!(c, m)
      @test c === rand!(rng, c, m)
      @test c === rand!(c, args...)
      @test c === rand!(rng, c, args...)
      @test reproducible(m)
      @test reproducible(args...)
    end
  end

  @testset "NumField/Coordinates" begin
    Qx, x = PolynomialRing(FlintQQ, "x")
    K, a = number_field(x^2+1, cached = false)
    BK = basis(K)
    for i = 1:5
      el = rand(K, -5:5)
      v = coordinates(el)
      @test el == dot(v, BK)
    end
    Kns, gKns = number_field([x^2+1, x^2+2], check = false)
    BKns = basis(Kns)
    for i = 1:5
      el = rand(Kns, -5:5)
      v = coordinates(el)
      @test el == dot(v, BKns)
    end
    Kt, t = PolynomialRing(K, "t", cached = false)
    L, gL = number_field(t^2+3, cached = false)
    BL = basis(L)
    BLabs = absolute_basis(L)
    for i = 1:5
      el = rand(L, -5:5)
      v = coordinates(el)
      vabs = absolute_coordinates(el)
      @test el == dot(v, BL)
      @test el == dot(vabs, BLabs)
    end
    Lns, gLns = number_field([t^2+3, t^2+5], check = false, cached = false)
    BLns = basis(Lns)
    BLnsabs = absolute_basis(Lns)
    for i = 1:5
      el = rand(Lns, -5:5)
      v = coordinates(el)
      vabs = absolute_coordinates(el)
      @test el == dot(v, BLns)
      @test el == dot(vabs, BLnsabs)
    end
    Knsy, y = PolynomialRing(Kns, "y", cached = false)
    F, gF = number_field(y^2+7, cached = false)
    BF = basis(F)
    BFabs = absolute_basis(F)
    for i = 1:5
      el = rand(F, -5:5)
      v = coordinates(el)
      vabs = absolute_coordinates(el)
      @test el == dot(v, BF)
      @test el == dot(vabs, BFabs)
    end
    Fns, gFns = number_field([y^2+7, y^2+11], cached = false)
    BFns = basis(Fns)
    BFnsabs = absolute_basis(Fns)
    for i = 1:5
      el = rand(Fns, -5:5)
      v = coordinates(el)
      vabs = absolute_coordinates(el)
      @test el == dot(v, BFns)
      @test el == dot(vabs, BFnsabs)
    end
    Fnsz, z = PolynomialRing(Fns, "z", cached = false)
    N, gN = number_field([z^2+17, z^2+19], check = false, cached = false)
    BN = basis(N)
    BNabs = absolute_basis(N)
    for i = 1:5
      el = rand(N, -5:5)
      v = coordinates(el)
      vabs = absolute_coordinates(el)
      @test el == dot(v, BN)
      @test el == dot(vabs, BNabs)
    end
  end

  @testset "relative extension" begin
    Qx, x = FlintQQ["x"]
    f = x^2 + 12x - 92
    K, a = NumberField(f, "a")
    Ky, y = K["y"]
    L, b = NumberField(y^2 + y + 1, "b")
    Lt, t = PolynomialRing(L)
    L1, gL1 = number_field([t^3-2])
    L1rel, mL1rel = relative_simple_extension(L1, K)
    el = mL1rel(gen(L1rel))
    @test isprimitive_over(el, K)

  end
end

@testset "issquare" begin
  #Enable one issquaer is working for non-monic defining polynomials
  Qx, x = QQ["x"]
  rangee = deleteat!(collect(-10:10), 11) # remove 0
  for d in [2,3,4,6,8,10]
    f = rand(Qx, d:d, -10:10)
    while !isirreducible(f)
      f = rand(Qx, d:d, -10:10)
    end
    K, a = NumberField(f)
    for m in 1:100
      b = rand(K, -10:10)//rand(rangee)
      c = b^2
      fl, d = issquare(c)
      @test fl
      @test d^2 == c
      b = rand(K, -10:10)
      KK, m = simplify(K)
      @test issquare(b)[1] == issquare(m\b)[1]
    end
  end

  f = -6//7*x + 1//9
  K, a = number_field(f, cached = false)
  c = 49//4
  fl, b = issquare(K(c))
  @test fl
  @test b^2 == c
end

@testset "Torsion units" begin
  Qx, x = QQ["x"]
  f = x^2 + 3
  K, a = NumberField(f, "a")
  G, mG = torsion_unit_group(K)
  @test order(G) == 6
  g = mG(G[1])
  @test g == torsion_units_generator(K)
  @test isone(g^6)
  @test !isone(g^3)
  @test !isone(g^2)
  @test Hecke.istorsion_unit_group_known(K)
  M = EquationOrder(K)
  A, mA = @inferred torsion_unit_group(M)
  @test order(A) == 2
  @test mA(A[1]) == M(-1)
  OK = maximal_order(K)
  A, mA = @inferred torsion_unit_group(OK)
  @test order(A) == 6
  @test isone(mA(A[1])^6)
  @test !isone(mA(A[1])^3)
  @test !isone(mA(A[1])^2)

  Kt, t = PolynomialRing(K, cached = false)
  Ls, gLs = number_field(t^2+1)
  G, mG = torsion_unit_group(Ls)
  @test ngens(G) == 1
  @test order(G) == 12
  @test torsion_units_order(Ls) == 12
  g = mG(G[1])
  @test g == torsion_units_generator(Ls)
  @test isone(g^12)
  @test !isone(g^4)
  @test !isone(g^3)

  Lns, gLns = number_field([t^2+1, t^2+2])
  G, mG = torsion_unit_group(Lns)
  @test ngens(G) == 1
  @test order(G) == 24
  @test torsion_units_order(Lns) == 24
  g = mG(G[1])
  @test g == torsion_units_generator(Lns)
  @test isone(g^24)
  @test !isone(g^8)
  @test !isone(g^3)
end


text:
using JuMP
import Ipopt

function P1(z)
    model = Model(Ipopt.Optimizer)
    @variable(model, x1)
    @variable(model, x2)
    @objective(model, Min, z[1] * x1 + z[2] * x2)
    @NLconstraint(model, c1, (x1 - 1)^2 + (x2 - 1)^2 <= 1)
    optimize!(model)
    return value.([x1, x2]), objective_value(model)
end
function P2(v)
    model = Model(Ipopt.Optimizer)
    @variable(model, z)
    @variable(model, x1)
    @variable(model, x2)
    @objective(model, Min, z)
    @NLconstraint(model, c1, (x1 - 1)^2 + (x2 - 1)^2 <= 1)
    @constraint(model, c2, x1 - z - v[1] <= 0)
    @constraint(model, c3, x2 - z - v[2] <= 0)
    optimize!(model)
    return (value.([x1, x2]), value(z))
end
function D2(v)
    model = Model(Ipopt.Optimizer)
    @variable(model, u >= 0)
    @variable(model, w1 >= 0)
    @variable(model, w2 >= 0)
    @NLobjective(
        model,
        Max,
        (-(w1^2 + w2^2) / (4 * u) - u) + w1 * (1 - v[1]) + w2 * (1 - v[2])
    )
    @constraint(model, c1, w1 + w2 == 1)
    optimize!(model)
    return (value.([w1, w2]), value(u))
end
function Γ(x)
    x
end


text:
using Printf, JuMP, DelimitedFiles, Mosek, MosekTools, LinearAlgebra

function read_data(INSTANCE_NAME, FORMULATION, DATA_PATH, REPO_PB)
  if FORMULATION == "NO_BLOCKS"
    BLOCKS=readdlm(joinpath(DATA_PATH,"decompositions","blocks_cholesky", INSTANCE_NAME*"_sdp_blocks.txt"))
    VAR = BLOCKS[:,2]
    B1 = ["B1" for i=1:size(BLOCKS,1)]
    BLOCKS = [B1 VAR]
    # CLIQUE_TREE = readdlm(joinpath(DATA_PATH, "cliquetree_MD",INSTANCE_NAME*"_sdp_cliquetree.txt"))
    CLIQUE_TREE = []
  else
    BLOCKS=readdlm(joinpath(DATA_PATH,"decompositions", "blocks_$FORMULATION",INSTANCE_NAME*"_sdp_blocks.txt"))
    CLIQUE_TREE = readdlm(joinpath(DATA_PATH,"decompositions", "cliquetree_$FORMULATION",INSTANCE_NAME*"_sdp_cliquetree.txt"))
  end
  CONSTRAINTS = readdlm(joinpath(DATA_PATH, REPO_PB, INSTANCE_NAME*"_sdp_ctr.txt"),skipstart=1)
  MATRICES = readdlm(joinpath(DATA_PATH, REPO_PB, INSTANCE_NAME*"_sdp_mat.txt"),skipstart=1)
  return BLOCKS, CLIQUE_TREE, CONSTRAINTS, MATRICES
end

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

function treat_constants(MATRICES_with_CONSTANTS, CONSTRAINTS)
  #find constant rows
  index_none_column2 = findall(x -> x == "NONE", MATRICES_with_CONSTANTS[:,2])
  index_none_column3 = findall(x -> x == "NONE", MATRICES_with_CONSTANTS[:,3])
  rows_to_remove = intersect(index_none_column2,index_none_column3)
  CONSTANTS = MATRICES_with_CONSTANTS[rows_to_remove,:]
  #remove constants from MATRICES
  MATRICES_without_CONSTANTS = MATRICES_with_CONSTANTS[setdiff(1:end, rows_to_remove), :]
  #define dict constant per objctr
  constants_per_objconstraint = Dict(CONSTRAINTS[i,1] => 0.0 for i=1:size(CONSTRAINTS, 1))
  constants_per_objconstraint["OBJ"] = 0.0
  for l=1:size(CONSTANTS,1)
    name = CONSTANTS[l,1]
    constants_per_objconstraint[name] = CONSTANTS[l,4]
  end
  return MATRICES_without_CONSTANTS, constants_per_objconstraint
end #function treat_constants

#########################################################################################################################
# liste des variables de chaque bloc
function create_list_of_variables_per_block(BLOCKS)
  DIFFERENT_BLOCKS = Set{String}(BLOCKS[:,1])
  block_var = Dict{String, Set{String}}()
  for block in DIFFERENT_BLOCKS
    block_var[block] = Set{String}()
  end
  for line=1:size(BLOCKS,1)
    block = BLOCKS[line,1]
    var = BLOCKS[line,2]
    block_var[block] = union(block_var[block],Set{String}([var]))
  end
  return block_var
end

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

function creating_and_solving_model(DIFFERENT_VAR, NB_BLOCKS, DIFFERENT_BLOCKS, block_var,CLIQUE_TREE, CONSTRAINTS, MATRICES, constants_per_constraint,my_solver, objscale)
  #initialize
  m = Model(with_optimizer(my_solver))
  NB_CONSTRAINTS = size(CONSTRAINTS,1)
  jumpX = Array{LinearAlgebra.Symmetric{VariableRef,Array{VariableRef,2}}}(undef,NB_BLOCKS)
  coeff_block = Dict{Tuple{String,String}, Set{String}}()
  mat_var = Dict{Tuple{String,String}, Dict{String,Any}}()
  println("Dictionary linking blocks with their set of variables created")
  for (block, var_list) in block_var
    for var1 in var_list
      for var2 in var_list
        mat_var[(var1,var2)] = Dict{String,Any}()
        coeff_block[(var1,var2)] = Set{String}()
      end
    end
  end
  println("Dictionary linking coeff with their couple of variables created")
  @printf("mat_var size is %8d\n", length(coeff_block))
  size_block = Dict{String, Int64}()
  i_block = 0
  for (block, var_list) in block_var
    i_block += 1
    i_var1 = 1
    size_block[block] = length(var_list)
    jumpX[i_block] =  @variable(m,[1:size_block[block],1:size_block[block]], base_name= "X_$block", PSD)
    for var1 in var_list
      i_var2 = 1
      for var2 in var_list
        mat_var[(var1,var2)][block] = jumpX[i_block][i_var1, i_var2]
        coeff_block[(var1,var2)] = union(coeff_block[(var1,var2)],Set{String}([block]))
        i_var2+=1
      end
      i_var1+=1
    end
  end
  println("Block variables defined in the model m : ", i_block, " - ", length(block_var))
  println("Dictionary linking variables, blocks and data coefficients done")
  xp = Dict{String,JuMP.GenericAffExpr}()
  xp["OBJ"] = 0*jumpX[1][1,1]
  for p = 1:NB_CONSTRAINTS
    xp[CONSTRAINTS[p,1]] = 0*jumpX[1][1,1]
  end
  NB_TERMS = size(MATRICES,1)
   #NO MEAN version
   for term = 1:NB_TERMS
     objctr = MATRICES[term,1]
     var1 = MATRICES[term,2]
     var2 = MATRICES[term,3]
     val = MATRICES[term,4]
     nb_block = length(coeff_block[(var1,var2)])
     min_size = Inf
     min_var = first(mat_var[(var1,var2)])[2]
     for (block, var) in mat_var[(var1,var2)]
       if length(block_var[block]) < min_size
         min_size = length(block_var[block])
         min_var = var
       end
     end
     add_to_expression!(xp[objctr], val*min_var)
   end
  println("Objective and constraints scalar products computed")
  my_timer = @elapsed @objective(m, Min , objscale*(xp["OBJ"]+constants_per_constraint["OBJ"]))
  @printf("%-35s%10.6f s\n", "Objective added to model m", my_timer)
  constraints_ref=Dict{String,JuMP.ConstraintRef}()
  for p = 1:NB_CONSTRAINTS
    name = CONSTRAINTS[p,1]
    if (CONSTRAINTS[p,2]=="EQ")
      constraints_ref[name] = @constraint(m,xp[name]+constants_per_constraint[name] == CONSTRAINTS[p,3])
    elseif (CONSTRAINTS[p,2]=="LEQ")
      constraints_ref[name] = @constraint(m,xp[name]+constants_per_constraint[name]<=CONSTRAINTS[p,4])
    elseif (CONSTRAINTS[p,2]=="GEQ")
      constraints_ref[name] = @constraint(m,xp[name]+constants_per_constraint[name]>=CONSTRAINTS[p,3])
    else
      constraints_ref[*(name,"_LEQ")] = @constraint(m,xp[name]+constants_per_constraint[name]<=CONSTRAINTS[p,4])
      constraints_ref[*(name,"_GEQ")] = @constraint(m,CONSTRAINTS[p,3]<=xp[name]+constants_per_constraint[name])
    end
  end
  println("initial constraints added to model m")
  #constraints linking common terms in blocks
  nb_coupling_constraints=0
  if NB_BLOCKS > 1
    for i in 1:size(CLIQUE_TREE,1)
      B1 = CLIQUE_TREE[i,1]
      B2 = CLIQUE_TREE[i,2]
      vars_B1 = block_var[B1]
      vars_B2 = block_var[B2]
      common_vars = [ var for var in intersect(vars_B1,vars_B2)]
      for i in 1:length(common_vars)
        var1 = common_vars[i]
        for j in i:length(common_vars)
          var2 = common_vars[j]
          JuMPvar1 = mat_var[(var1,var2)][B1]
          JuMPvar2 = mat_var[(var1,var2)][B2]
          constraints_ref[*("cc$(B1)_$(B2)_$var1","_$var2")]=@constraint(m,JuMPvar1-JuMPvar2==0)
          nb_coupling_constraints+=1
        end
      end
    end
  end
  println("nb_coupling_constraints primal : ", nb_coupling_constraints)
  println("nb variables primal : ", length(JuMP.all_variables(m)))
  for ctr_type in list_of_constraint_types(m)
    println("nb constraints type $ctr_type primal : ", num_constraints(m, ctr_type[1], ctr_type[2]))
  end
  # print(m)
  val, time, bytes, gctime, memallocs = @timed optimize!(m)
  # MOI.write_to_file(backend(m).optimizer.model, "dump.task.gz")
  mem = bytes*10^(-6)
  println("Resolution time : ", time)
  println("Memory consumption $mem Mo")
  println("\nObjective value: ", 1/objscale*JuMP.objective_value(m))
  Xoptimal=Array{Any}(nothing,NB_BLOCKS)
  for i=1:NB_BLOCKS
    Xoptimal[i]=JuMP.value.(jumpX[i])
  end
  return Xoptimal,JuMP.objective_value(m), JuMP.termination_status(m), time, mem
end

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

function sdp_model_solver_to_specify(INSTANCE_NAME, FORMULATION, DATA_PATH, REPO, my_solver, objscale)
  my_solver = Mosek.Optimizer
  #read data files
  my_timer=@elapsed (BLOCKS, CLIQUE_TREE, CONSTRAINTS, MATRICES_with_CONSTANTS)=read_data(INSTANCE_NAME, FORMULATION, DATA_PATH, REPO)
  @printf("%-35s%10.6f s\n", "read_data", my_timer)
  #treat constants
  my_timer = @elapsed (MATRICES, constants_per_objconstraint) = treat_constants(MATRICES_with_CONSTANTS, CONSTRAINTS)
  @printf("%-35s%10.6f s\n", "constant coeffcients treated", my_timer)

  #CONSTRAINTS
  DIFFERENT_CONSTRAINTS = Set{String}(CONSTRAINTS[:,1])
  NB_CONSTRAINTS = length(DIFFERENT_CONSTRAINTS)
  #variables
  DIFFERENT_VAR = Set{String}(BLOCKS[:,2])
  NB_VAR = length(DIFFERENT_VAR)
  #define block patterns
  DIFFERENT_BLOCKS = Set{String}(BLOCKS[:,1])
  NB_BLOCKS = length(DIFFERENT_BLOCKS)
  println("\n NB BLOCKS = $NB_BLOCKS \n")
  my_timer = @elapsed block_var = create_list_of_variables_per_block(BLOCKS)
  @printf("%-35s%10.6f s\n", "create_list_of_variables_per_block", my_timer)
  my_timer= @elapsed (Xoptimal,obj,status, time, mem) = creating_and_solving_model(DIFFERENT_VAR, NB_BLOCKS, DIFFERENT_BLOCKS, block_var,CLIQUE_TREE, CONSTRAINTS, MATRICES, constants_per_objconstraint,my_solver, objscale)
  #return S, Xoptimal, y
  @printf("%-35s%10.6f s\n", "creating_and_solving_model", my_timer)
  return Xoptimal, obj, status, time, mem
end

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

function solve_sdp(INSTANCE_NAME, FORMULATION, DATA_PATH)
    objscale = 4
    println("objscale : 10^(-$objscale) \n")
    Xoptimal, obj, status, time, mem = sdp_model_solver_to_specify(INSTANCE_NAME, FORMULATION, DATA_PATH,"OPF_instances", Mosek.Optimizer, 10.0^(-objscale))
    return
end

max_stars_repo_name: JulieSliwak/CliqueMerging

max_issues_repo_name: JulieSliwak/CliqueMerging

max_forks_repo_name: JulieSliwak/CliqueMerging


text:
module InvPendulum

    include("src/invpend_julia.jl")
    rk4_func  = inv_pendulum_rk4  # set standard Environment function
    cost_func = inv_pendulum_cost#_swingup # set standard Environment cost-function
    s_init = [0.0, 0.0, 0.0, 0.0] # set standard Environment initial position
    #s_init = [0.0, 0.0, 0.0, pi]
    u_range = [11.1]
    export rk4_func, cost_func, s_init, u_range

end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
isempty(@__DIR__) || push!(LOAD_PATH, joinpath(@__DIR__, "..", "src"))
using Render

function main()
    if length(ARGS) < 1
        include(joinpath(@__DIR__, "repl.jl"))
    else
        map(x -> open(Render.renderfile, x, "r"), ARGS)
    end
end

isinteractive() || @time main()

max_stars_repo_name: AnAverageHuman/jlgraphics

max_issues_repo_name: AnAverageHuman/jlgraphics

max_forks_repo_name: AnAverageHuman/jlgraphics


text:
macro eshead(info,  url , query  )

	esc(
		quote
			try
				respos = HTTP.request("HEAD", HTTP.URI($(url)), $(info).jheader , query= $query ;  $(info).conf...)

				if respos.status == 200
					"OK"
				end
			catch
				"Fail"
			end

		end )
end

macro esdelete(info,  url,   query   )
	esc(
		quote
			respos = HTTP.request("DELETE", HTTP.URI($(url)) , $(info).jheader,  query= $query ;  $(info).conf...)

			if respos.status == 200
				JSON.parse(String(respos.body))
			end

		end )
end
macro catexport(info, method, url , query  )
	esc(
		quote
			respos = HTTP.request($method, HTTP.URI($(url)), $(info).jheader, query= $query ;  $(info).conf...)
			String(respos.body) |> println
		end )
end

macro esexp(x,y)
	esc(:(get($x, $y, missing)))
end

macro esexport(info, method, url, body , query , header )
	esc(
		quote

			respos = HTTP.request($method, HTTP.URI($(url)) , $header , $body, query= $query;  $(info).conf...)

			if respos.status == 200
				JSON.parse(String(respos.body))
			end

		end )
end

function genfunction(  kw::Vector  )

	func = Expr(:call ,
		kw[2] ,
		Expr(:parameters,
			Expr(:(...), :kw)),
		Expr(:(::) , :info , :Esinfo)
		)

	url = :( url = make_url( $(kw[3]), info ) )

	if length(kw) >= 5
		append!( url.args[2].args ,  kw[5:end] )
		append!(func.args, Expr.(:(::) , kw[5:end] , :AbstractString))
	end

	if kw[4] >= 1
		push!(func.args,  Expr(:(::) , :body , :T) )
		body = :body
	else
		body = Dict()
	end

	block = Expr(:block ,
		:( querys = Dict(kw...) ),
		url ,
		Expr(:macrocall ,
				 Symbol("@esexport") ,
				 "",
				 :info,
				 kw[1],
				 :url ,
				 Expr(:call , :ifelse , Expr(:call, :isa, body, :Dict),
							Expr(:call, :json, body) , body ) ,
				 :querys ,
				 Expr(:(.), :info , :(:jheader))
				 )
			)
	if kw[4] >= 1
		esc(Expr(:function, Expr(:where , func,
					Expr(:(<:) , :T , Expr(:curly,:Union, :Dict, :AbstractString))), block))
	else
		esc(Expr(:function, func, block))
	end

end

macro genfunction( kw... )

	 genfunction( collect(kw) )

end

max_stars_repo_name: frankill/ES

max_issues_repo_name: frankill/ES

max_forks_repo_name: frankill/ES


text:
using MyTestPackage
using Test

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

max_stars_repo_name: blixi/MyTestPackage

max_issues_repo_name: blixi/MyTestPackage

max_forks_repo_name: blixi/MyTestPackage


text:
using Flux
using Flux: @epochs
using Statistics
using Random
using Parameters: @with_kw
using IterTools: ncycle

@with_kw mutable struct Args
    lr::Float64 = 0.5
    batch_size::Int = 256
    repeat::Int = 20
end

images = Flux.Data.FashionMNIST.images()
labels = Flux.Data.FashionMNIST.labels()

function get_fashion_label(labels)
	text_label = ["t-shirt", "trouser", "pullover", "dress", "coat", "sandal", "shirt", "sneaker", "bag", "ankle boot"]
	return [text_label[i+1] for i in labels]
end

#get train_data

data_x = rand(784, 6000)
data_y = []

for i in 1:6000
	data_x[:,i] = Float64.(reshape(images[i],(784,1)))
	push!(data_y, get_fashion_label(labels)[i])
end

label = sort(unique(data_y))
data_onehot_labels = Flux.onehotbatch(data_y, label)

train_x = data_x[:, [1:3:6000 ; 2:3:6000]]
train_y = data_onehot_labels[:, [1:3:6000 ; 2:3:6000]]

test_x = data_x[:, 3:3:6000]
test_y = data_onehot_labels[:, 3:3:6000]

train_data = Flux.Data.DataLoader((train_x, train_y), batchsize=Args().batch_size, shuffle=true)

#model
model = Chain(
    Dense(784, 256, relu),
    Dense(256, 10)
    )

#define loss function: cross entropy function
loss(x, y) = Flux.logitcrossentropy(model(x), y)

# params
ps = Flux.params(model)

#SDG
opt = Descent(Args().lr)

@time Flux.train!(loss, ps, ncycle(train_data, Args().repeat),opt)

accuracy(x, y, model) = Flux.mean(Flux.onecold(model(x)) .== Flux.onecold(y))

print("train loss: ", accuracy(train_x, train_y, model), ", test loss:", accuracy(test_x, test_y, model))


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Integer API #
# ------------------------------------------------------------
initialize_integer_algorithm!(lshaped::AbstractLShaped) = initialize_integer_algorithm!(lshaped.integer, lshaped.structure.first_stage)
initialize_integer_algorithm!(subproblem::SubProblem) = initialize_integer_algorithm!(subproblem.integer_algorithm, subproblem)
# ------------------------------------------------------------
# Attributes #
# ------------------------------------------------------------
"""
    RawIntegerAlgorithmParameter

An optimizer attribute used for raw parameters of the integer algorithm. Defers to `RawParameter`.
"""
struct RawIntegerParameter <: IntegerParameter
    name::Any
end

include("common.jl")
include("ignore_integers.jl")
include("combinatorial_cuts.jl")
include("convexification.jl")


text:
using Yota
import Yota: TReal, TAny, record!, Call, grad!
using XGrad
using BenchmarkTools

include("functions.jl")

# xgrad derivatives
@diffrule logistic(x::Number) x (logistic(x) .* (1 .- logistic(x)) .* ds)

# yota wrappers and derivatives
logistic(x::TReal) = record!(x.tape, Call, logistic, (x,))

function grad!(dy::TAny, ::Val{1}, op::Call{typeof(logistic), Tuple{TReal}})
    x = op.args[1]
    ll = logistic(x)
    return (ll * (1 - ll) * dy)
end

# runner

function perf_test(f; ctx=Dict(), compile_tape=true, inputs...)
    vals = ([val for (name, val) in inputs]...,)

    println("Compiling derivatives using XGrad")
    @time df = xdiff(f; ctx=ctx, inputs...)
    mem = Dict()
    println("Testing XGrad...")
    r1 = @benchmark $df($vals...; mem=$mem)
    show(stdout, MIME{Symbol("text/plain")}(), r1)
    println("\n")

    println("Compiling derivatives using Yota")
    @time grad(f, vals...)
    r2 = @benchmark grad($f, $(vals)...)

    show(stdout, MIME{Symbol("text/plain")}(), r2)
    println("\n")

end

# benchmarks

function benchmark_autoencoder()
    f = autoencoder_cost
    println("\n## On larger data\n")
    We1 = rand(2000, 10_000); b1 = rand(2000); We2 = rand(1000, 2000); b2 = rand(1000);
    Wd = rand(10_000, 1000); x = rand(10_000, 100);
    inputs = [:We1 => We1, :We2 => We2, :Wd => Wd, :b1 => b1, :b2 => b2, :x => x];
    perf_test(f; inputs...)

    println("\n## On smaller data\n")
    We1 = rand(200, 1000); b1 = rand(200); We2 = rand(100, 200); b2 = rand(100);
    Wd = rand(1000, 100); x = rand(1000, 100);
    inputs = [:We1 => We1, :We2 => We2, :Wd => Wd, :b1 => b1, :b2 => b2, :x => x];
    perf_test(f; compile_tape=false, inputs...)
end

function benchmark_mlp1()
    f = mlp1
    println("\n## On larger data\n")
    w1=rand(2000, 10000); w2=rand(1000, 2000); w3=rand(1000, 1000); x1=rand(10000, 500);
    inputs = [:w1=>w1, :w2=>w2, :w3=>w3, :x1=>x1];
    perf_test(f; inputs...)

    println("\n## On smaller data\n")
    w1=rand(200, 1000); w2=rand(100, 200); w3=rand(100, 100); x1=rand(1000, 10);
    inputs = [:w1=>w1, :w2=>w2, :w3=>w3, :x1=>x1];
    perf_test(f; inputs...)
end


text:
module Audios

using FileType
TypeMidi = FileType.Types("mid", MIME("audio/midi"))
TypeMp3  = FileType.Types("mp3", MIME("audio/mpeg"))
TypeM4a  = FileType.Types("m4a", MIME("audio/m4a"))
TypeOgg  = FileType.Types("ogg", MIME("audio/ogg"))
TypeFlac = FileType.Types("flac", MIME("audio/x-flac"))
TypeWav  = FileType.Types("wav", MIME("audio/x-wav"))
TypeAmr  = FileType.Types("amr", MIME("audio/amr"))
TypeAac  = FileType.Types("aac", MIME("audio/aac"))

function Midi(buf)::Bool
	return length(buf) > 3 &&
		buf[1] == 0x4D && buf[2] == 0x54 &&
		buf[3] == 0x68 && buf[4] == 0x64
end

function Mp3(buf)::Bool
	return length(buf) > 2 &&
		((buf[1] == 0x49 && buf[2] == 0x44 && buf[3] == 0x33) ||
			(buf[1] == 0xFF && buf[2] == 0xfb))
end

function M4a(buf)::Bool
	return length(buf) > 10 &&
		((buf[5] == 0x66 && buf[6] == 0x74 && buf[7] == 0x79 &&
			buf[8] == 0x70 && buf[9] == 0x4D && buf[10] == 0x34 && buf[11] == 0x41) ||
			(buf[1] == 0x4D && buf[2] == 0x34 && buf[3] == 0x41 && buf[4] == 0x20))
end

function Ogg(buf)::Bool
	return length(buf) > 3 &&
		buf[1] == 0x4F && buf[2] == 0x67 &&
		buf[3] == 0x67 && buf[4] == 0x53
end

function Flac(buf)::Bool
	return length(buf) > 3 &&
		buf[1] == 0x66 && buf[2] == 0x4C &&
		buf[3] == 0x61 && buf[4] == 0x43
end

function Wav(buf)::Bool
	return length(buf) > 11 &&
		buf[1] == 0x52 && buf[2] == 0x49 &&
		buf[3] == 0x46 && buf[4] == 0x46 &&
		buf[9] == 0x57 && buf[10] == 0x41 &&
		buf[11] == 0x56 && buf[12] == 0x45
end

function Amr(buf)::Bool
	return length(buf) > 11 &&
		buf[1] == 0x23 && buf[2] == 0x21 &&
		buf[3] == 0x41 && buf[4] == 0x4D &&
		buf[5] == 0x52 && buf[6] == 0x0A
end

function Aac(buf)::Bool
	return length(buf) > 1 &&
		((buf[1] == 0xFF && buf[2] == 0xF1) ||
			(buf[1] == 0xFF && buf[2] == 0xF9))
end

Audio = Dict(
	TypeMidi => Midi,
	TypeMp3 =>  Mp3,
	TypeM4a => M4a,
	TypeOgg =>  Ogg,
	TypeFlac => Flac,
	TypeWav =>  Wav,
	TypeAmr =>  Amr,
	TypeAac =>  Aac,
)

end


text:
function simulate(	model::String,
					m::Array{Float64, 1},
					c::Array{Float64, 1},
					ax::Float64,
					ay::Float64,
					az::Array{Float64, 1},
					Lx::Float64,
					Ly::Float64,
					Lz::Float64,
					number_of_frames::Array{Int64, 1},
					deltat::Float64,
					kmin::Int64,
					number_of_de_bins::Int64,
					ub_de::Float64)

	# Number of components in distribution.
	number_of_components::Int64 = length(m)
	fractions::Array{Float64, 1} = c / sum(c)
	cum_fractions::Array{Float64, 1} = cumsum(fractions)

	# Intensity of Poisson distribution of the number of particles. The factor
	# 1e12 takes into account that concentration is specified in particles/ml.
	lambda::Float64 = sum(c) * Lx * Ly * Lz / 1e12
	kmax::Int64 = maximum(number_of_frames)

	# Number of videos.
	number_of_videos::Int64 = length(number_of_frames)

	# Lower and upper bounds for the detection region.
	lbx::Float64 = 0.5 * (Lx - ax)
	ubx::Float64 = 0.5 * (Lx + ax)
	lby::Float64 = 0.5 * (Ly - ay)
	uby::Float64 = 0.5 * (Ly + ay)
	lbz::Array{Float64, 1} = 0.5 * (Lz - az)
	ubz::Array{Float64, 1} = 0.5 * (Lz + az)

	# Variables for storing positions, displacements, number of particles, etc
	x::Float64 = 0.0
	y::Float64 = 0.0
	z::Float64 = 0.0
	deltax::Float64 = 0.0
	deltay::Float64 = 0.0
	deltaz::Float64 = 0.0
	number_of_particles::Int64 = 0
	k::Int64 = 0
	de::Float64 = 0.0
	std_dev_random_walk::Float64 = 0.0
	index::Int64 = 0

	# Histogram matrix for number of positions and estimated diffusion coefficients of all
	# recorded trajectories.
	H::Array{Int64, 2} = zeros(kmax, number_of_de_bins)
	d_de::Float64 = ub_de / convert(Float64, number_of_de_bins)
	idx::Int64 = 0

	for current_video = 1:number_of_videos
		number_of_particles = rand_poisson(lambda)

		for current_particle = 1:number_of_particles
			# Generate random diffusion coefficent from distribution, or more precisely,
			# a random standard deviation for the displacements.
			index = rand_weighted_index(cum_fractions)
			std_dev_random_walk = sqrt(2.0 * m[index] * deltat)

			# Random initial position.
			x = Lx * rand()
			y = Ly * rand()
			z = Lz * rand()

			# Reset empirical diffusion coefficient.
			de = 0.0

			# In detection region or not? Check z limits first because they are more restrictive.
			#if (lbz[index] <= z <= ubz[index]) & (lbx <= x <= ubx) & (lby <= y <= uby)
			if (abs(z - 0.5 * Lz) <= 0.5 * az[index]) & (abs(x - 0.5 * Lx) <= 0.5 * ax) & (abs(y - 0.5 * Ly) <= 0.5 * ay)
				k = 1
			else
				k = 0
			end

			# Let particle diffuse through remainder of video.
			for current_frame = 2:number_of_frames[current_video]
				deltax = std_dev_random_walk * randn()
				deltay = std_dev_random_walk * randn()
				deltaz = std_dev_random_walk * randn()

				x = x + deltax
				y = y + deltay
				z = z + deltaz

				x = position_periodic(x, Lx)
				y = position_periodic(y, Ly)
				z = position_periodic(z, Lz)

				#if (lbz[index] <= z <= ubz[index]) & (lbx <= x <= ubx) & (lby <= y <= uby)
				if (abs(z - 0.5 * Lz) <= 0.5 * az[index]) & (abs(x - 0.5 * Lx) <= 0.5 * ax) & (abs(y - 0.5 * Ly) <= 0.5 * ay)
					k = k + 1
					if k >= 2 # The first point inside the detection region does not contributed toward the estimated diffusion coefficient.
						de = de + deltax^2 + deltay^2 # Only in x-y plane.
					end
				elseif k > 0
					if k >= kmin
						de = de / (convert(Float64, k - 1) * 4.0 * deltat) # The '4' comes from the 2D observations.
						idx = convert(Int64, ceil(de / d_de))
						if 1 <= idx <= number_of_de_bins
							H[k, idx] += 1
						end
					end
					k = 0
					de = 0.0
				end
			end

			# If particle is inside detection region in last frame, we have now missed that trajectory and need to add it now.
			if k >= kmin # Assuming kmin >= 2 otherwise the division fails.
				de = de / (convert(Float64, k - 1) * 4.0 * deltat) # The '4' comes from the 2D observations.
				idx = convert(Int64, ceil(de / d_de))
				if 1 <= idx <= number_of_de_bins
					H[k, idx] += 1
				end
			end

		end
	end

	return H
end


text:
"""
`SaveModel` module purpose is to provide a delegate-to-solver option to
`MathOptInterface.write_to_file`. Instead of using the current MOI
fallback (that only works for MPS and is kinda memory-hungry and slow),
`SaveModel.write_to_file` checks if `SaveModel.natively_supports` returns
`true` for its arguments, and in such case calls
`SaveModel.native_write_to_file` over the object returned by the
`MOI.RawSolver()` property; otherwise it calls the `MOI` fallback.

"""
module SaveModel

import JuMP
import MathOptInterface
const MOI = MathOptInterface

import Requires.@require

import TimerOutputs.@timeit
import ..TIMER

# Based on: https://discourse.julialang.org/t/trying-to-save-mps-from-raw-solver-and-problems-with-moi-rawsolver/58515/2?u=henrique_becker
_get_raw_solver(model :: JuMP.Model) = _get_raw_solver(JuMP.backend(model))
function _get_raw_solver(model :: MOI.Utilities.CachingOptimizer)
	if model.state == MOI.Utilities.EMPTY_OPTIMIZER
		MOI.Utilities.attach_optimizer(model)
	end
	return _get_raw_solver(model.optimizer)
end
function _get_raw_solver(model :: MOI.Bridges.LazyBridgeOptimizer)
	_get_raw_solver(model.model)
end
_get_raw_solver(model) = model # Unrecognized means we arrived were we wanted.

"""
	write_to_file(model, filename; format = FORMAT_AUTOMATIC

Save the model to filename trying to use the solver mechanism instead of
the generic `JuMP`/`MathOptInterface` one.
"""
@timeit TIMER function write_to_file(
	model :: JuMP.Model,
	filename :: String;
	format :: MOI.FileFormats.FileFormat = MOI.FileFormats.FORMAT_AUTOMATIC
)
	raw_solver = _get_raw_solver(model)
	if natively_supports(raw_solver, filename; format = format)
		native_write_to_file(raw_solver, filename; format = format)
	else
		JuMP.write_to_file(model, filename; format = format)
	end

	return
end
export write_to_file

function natively_supports(
	raw_solver, filename :: String;
	format :: MOI.FileFormats.FileFormat = MOI.FileFormats.FORMAT_AUTOMATIC
) :: Bool
	return false
end

function native_write_to_file(
	raw_solver, filename :: String;
	format :: MOI.FileFormats.FileFormat = MOI.FileFormats.FORMAT_AUTOMATIC
)
	throw(ArgumentError(
		"`native_write_to_file` was called for a solver it does know." *
		" Check if the the solver package is loaded, and if some package" *
		" provides a solver-specific method for the functions in" *
		" `GuillotineModels.SaveModel`."
	))
end

const MOI_FF = MOI.FileFormats
const MOI_FF_FF = MOI_FF.FileFormat

const _KNOWN_FORMATS = Tuple{String, MOI_FF_FF}[
	(".cbf", MOI_FF.FORMAT_CBF),
	(".lp", MOI_FF.FORMAT_LP),
	(".mof.json", MOI_FF.FORMAT_MOF),
	(".mps", MOI_FF.FORMAT_MPS),
	(".sdpa", MOI_FF.FORMAT_SDPA),
	(".dat-s", MOI_FF.FORMAT_SDPA),
]

function _format_from_ext(filename :: String)
	for (ext, enum) in _KNOWN_FORMATS
		(endswith(filename, ext) || occursin("$ext.", filename)) && return enum
	end
	error("Unable to automatically detect format of $(filename).")
end

function _real_format(filename :: String, format :: MOI_FF_FF) :: MOI_FF_FF
	return if format == MOI_FF.FORMAT_AUTOMATIC
		_format_from_ext(filename) :: MOI_FF_FF
	else
		format
	end
end

function __init__()
	@require CPLEX="a076750e-1247-5638-91d2-ce28b192dca0" begin
		function natively_supports(
			model :: CPLEX.Optimizer, filename :: String;
			format :: MOI_FF_FF = MOI_FF.FORMAT_AUTOMATIC
		)
			real_format = _real_format(filename, format)
			return real_format in (MOI_FF.FORMAT_MPS, MOI_FF.FORMAT_LP)
		end

		function native_write_to_file(
			model :: CPLEX.Optimizer, filename :: String;
			format :: MOI_FF_FF = MOI_FF.FORMAT_AUTOMATIC
		)
			@assert format == MOI_FF.FORMAT_AUTOMATIC ||
				format == _format_from_ext(filename)
			return CPLEX.CPXwriteprob(model.env, model.lp, filename, C_NULL)
		end
	end # @require CPLEX

	@require Gurobi="2e9cd046-0924-5485-92f1-d5272153d98b" begin
		function natively_supports(
			model :: Gurobi.Optimizer, filename :: String;
			format :: MOI_FF_FF = MOI_FF.FORMAT_AUTOMATIC
		)
			real_format = _real_format(filename, format)
			return real_format in (MOI_FF.FORMAT_MPS, MOI_FF.FORMAT_LP)
		end
		function native_write_to_file(
			model :: Gurobi.Optimizer, filename :: String;
			format :: MOI_FF_FF = MOI_FF.FORMAT_AUTOMATIC
		)
			@assert format == MOI_FF.FORMAT_AUTOMATIC ||
				format == _format_from_ext(filename)
			# Calling GRBupdatemodel directly here is a recipe for a disaster,
			# we need to use Gurobi.jl internals to avoid any subtle bugs.
			Gurobi._update_if_necessary(model)
			return Gurobi.GRBwrite(model, filename)
		end
	end # @require Gurobi
end # __init__()

end # module


text:
@testset "Label propagation" begin
    n = 10
    g10 = CompleteGraph(n)
    for g in testgraphs(g10)
        z = copy(g)
        for k = 2:5
            z = blockdiag(z, g)
            add_edge!(z, (k - 1) * n, k * n)
            c, ch = @inferred(label_propagation(z))
            a = collect(n:n:(k * n))
            a = Int[div(i - 1, n) + 1 for i = 1:(k * n)]
          # check the number of communities
            @test length(unique(a)) == length(unique(c))
          # check the partition
            @test a == c
        end
    end
end


text:
using Test
using MPI: mpiexec

# run tests on Travis CI in parallel
const TRIXI_TEST = get(ENV, "TRIXI_TEST", "all")
const ON_APPVEYOR = lowercase(get(ENV, "APPVEYOR", "false")) == "true"
const TRIXI_MPI_NPROCS = clamp(Sys.CPU_THREADS, 2, 3)

@time @testset "Trixi.jl tests" begin
  @time if TRIXI_TEST == "all" || TRIXI_TEST == "1D"
    include("test_examples_1d.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "2D"
    include("test_examples_2d.jl")
  end

  @time if (TRIXI_TEST == "all" && !ON_APPVEYOR) || TRIXI_TEST == "3D"
    include("test_examples_3d.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "misc"
    include("test_manual.jl")
    include("test_elixirs.jl")
  end

  @time if (TRIXI_TEST == "all" && !ON_APPVEYOR) || TRIXI_TEST == "paper-self-gravitating-gas-dynamics"
    include("test_paper-self-gravitating-gas-dynamics.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "parallel_2d"
    # Based on `runtests.jl` from `MPI.jl` and `PencilArrays.jl`
    # Precompilation disabled to prevent race conditions when loading packages
    mpiexec() do cmd
      run(`$cmd -n $TRIXI_MPI_NPROCS $(Base.julia_cmd()) --compiled-modules=no test_examples_parallel_2d.jl`)
    end
  end
end


text:
export AbstractMessage, Message, VariationalMessage
export getdata, is_clamped, is_initial, as_message
export multiply_messages

using Distributions
using Rocket

import Rocket: getrecent
import Base: *, +, ndims, precision, length, size, show

"""
    AbstractMessage

An abstract supertype for all concrete message types.

See also: [`Message`](@ref)
"""
abstract type AbstractMessage end

"""
    materialize!(message::AbstractMessage)

Materializes an abstract message and converts it to be of type `Message`.

See also: [`Message`](@ref)
"""
function materialize! end

"""
    Message{D} <: AbstractMessage

`Message` structure encodes a **Belief Propagation** message, which holds some `data` that usually a probability distribution, but can also be an arbitrary object.
Message acts as a proxy structure to `data` object and proxies most of the statistical functions, e.g. `mean`, `mode`, `cov` etc.

# Arguments
- `data::D`: message always holds some data object associated with it
- `is_clamped::Bool`, specifies if this message is clamped
- `is_initial::Bool`, specifies if this message is initial

# Example

```jldoctest
julia> distribution = Gamma(10.0, 2.0)
Gamma{Float64}(α=10.0, θ=2.0)

julia> message = Message(distribution, false, true)
Message(Gamma{Float64}(α=10.0, θ=2.0))

julia> mean(message)
20.0

julia> getdata(message)
Gamma{Float64}(α=10.0, θ=2.0)

julia> is_clamped(message)
false

julia> is_initial(message)
true

```

See also: [`AbstractMessage`](@ref), [`materialize!`](@ref)
"""
struct Message{D} <: AbstractMessage
    data       :: D
    is_clamped :: Bool
    is_initial :: Bool
end

getdata(message::Message)    = message.data
is_clamped(message::Message) = message.is_clamped
is_initial(message::Message) = message.is_initial

getdata(messages::NTuple{ N, <: Message }) where N = map(getdata, messages)

materialize!(message::Message) = message

Base.show(io::IO, message::Message) = print(io, string("Message(", getdata(message), ")"))

Base.:*(left::Message, right::Message) = multiply_messages(ProdAnalytical(), left, right)

function multiply_messages(prod_parametrisation, left::Message, right::Message)
    # We propagate clamped message, in case if both are clamped
    is_prod_clamped = is_clamped(left) && is_clamped(right)
    # We propagate initial message, in case if both are initial or left is initial and right is clameped or vice-versa
    is_prod_initial = !is_prod_clamped && (is_clamped_or_initial(left)) && (is_clamped_or_initial(right))

    return Message(prod(prod_parametrisation, getdata(left), getdata(right)), is_prod_clamped, is_prod_initial)
end

# Note: we need extra Base.Generator(as_message, messages) step here, because some of the messages might be VMP messages
# We want to cast it explicitly to a Message structure (which as_message does in case of VariationalMessage)
# We use with Base.Generator to reduce an amount of memory used by this procedure since Generator generates items lazily
prod_foldl_reduce(prod_constraint, form_constraint, ::FormConstraintCheckEach) = (messages) -> foldl((left, right) -> constrain_form(form_constraint, multiply_messages(prod_constraint, left, right)), Base.Generator(as_message, messages))
prod_foldl_reduce(prod_constraint, form_constraint, ::FormConstraintCheckLast) = (messages) -> constrain_form(form_constraint, foldl((left, right) -> multiply_messages(prod_constraint, left, right), Base.Generator(as_message, messages)))

prod_foldr_reduce(prod_constraint, form_constraint, ::FormConstraintCheckEach) = (messages) -> foldr((left, right) -> constrain_form(form_constraint, multiply_messages(prod_constraint, left, right)), Base.Generator(as_message, messages))
prod_foldr_reduce(prod_constraint, form_constraint, ::FormConstraintCheckLast) = (messages) -> constrain_form(form_constraint, foldr((left, right) -> multiply_messages(prod_constraint, left, right), Base.Generator(as_message, messages)))

# Base.:*(m1::Message, m2::Message) = multiply_messages(m1, m2)

Distributions.pdf(message::Message, x)    = Distributions.pdf(getdata(message), x)
Distributions.logpdf(message::Message, x) = Distributions.logpdf(getdata(message), x)

MacroHelpers.@proxy_methods Message getdata [
    Distributions.mean,
    Distributions.median,
    Distributions.mode,
    Distributions.shape,
    Distributions.scale,
    Distributions.rate,
    Distributions.var,
    Distributions.std,
    Distributions.cov,
    Distributions.invcov,
    Distributions.logdetcov,
    Distributions.entropy,
    Distributions.params,
    Base.precision,
    Base.length,
    Base.ndims,
    Base.size,
    mean_cov,
    mean_var,
    mean_invcov,
    mean_precision,
    weightedmean_cov,
    weightedmean_var,
    weightedmean_invcov,
    weightedmean_precision,
    probvec,
    weightedmean,
    inversemean,
    logmean,
    meanlogmean,
    mirroredlogmean,
    loggammamean
]

## Variational Message

mutable struct VariationalMessage{R, S, F} <: AbstractMessage
    messages   :: R
    marginals  :: S
    mappingFn  :: F
    cache      :: Union{Nothing, Message}
end

VariationalMessage(messages::R, marginals::S, mappingFn::F) where { R, S, F } = VariationalMessage(messages, marginals, mappingFn, nothing)

Base.show(io::IO, ::VariationalMessage) = print(io, string("VariationalMessage(:postponed)"))

getcache(vmessage::VariationalMessage)                    = vmessage.cache
setcache!(vmessage::VariationalMessage, message::Message) = vmessage.cache = message

compute_message(vmessage::VariationalMessage) = vmessage.mappingFn((vmessage.messages, getrecent(vmessage.marginals)))

function materialize!(vmessage::VariationalMessage)
    cache = getcache(vmessage)
    if cache !== nothing
        return cache
    end
    message = compute_message(vmessage)
    setcache!(vmessage, message)
    return message
end

## Utility functions

as_message(message::Message)             = message
as_message(vmessage::VariationalMessage) = materialize!(vmessage)

## Operators

# TODO
reduce_messages(messages) = mapreduce(as_message, (left, right) -> multiply_messages(ProdAnalytical(), left, right), messages)

## Message Mapping structure
## https://github.com/JuliaLang/julia/issues/42559
## Explanation: Julia cannot fully infer type of the lambda callback function in activate! method in node.jl file
## We create a lambda-like callable structure to improve type inference and make it more stable
## However it is not fully inferrable due to dynamic tags and variable constraints, but still better than just a raw lambda callback

struct MessageMapping{F, T, C, N, M, A, R}
    vtag            :: T
    vconstraint     :: C
    msgs_names      :: N
    marginals_names :: M
    meta            :: A
    factornode      :: R
end

message_mapping_fform(::MessageMapping{F}) where F = F
message_mapping_fform(::MessageMapping{F}) where F <: Function = F.instance

function MessageMapping(::Type{F}, vtag::T, vconstraint::C, msgs_names::N, marginals_names::M, meta::A, factornode::R) where { F, T, C, N, M, A, R }
    return MessageMapping{F, T, C, N, M, A, R}(vtag, vconstraint, msgs_names, marginals_names, meta, factornode)
end

function MessageMapping(::F, vtag::T, vconstraint::C, msgs_names::N, marginals_names::M, meta::A, factornode::R) where { F <: Function, T, C, N, M, A, R}
    return MessageMapping{F, T, C, N, M, A, R}(vtag, vconstraint, msgs_names, marginals_names, meta, factornode)
end

function (mapping::MessageMapping)(dependencies)

    messages  = dependencies[1]
    marginals = dependencies[2] # getrecent(marginals) call happens in VariationalMessage

    # Message is clamped if all of the inputs are clamped
    is_message_clamped = __check_all(is_clamped, messages) && __check_all(is_clamped, marginals)

    # Message is initial if it is not clamped and all of the inputs are either clamped or initial
    is_message_initial = !is_message_clamped && (__check_all(is_clamped_or_initial, messages) && __check_all(is_clamped_or_initial, marginals))

    message = rule(
        message_mapping_fform(mapping),
        mapping.vtag,
        mapping.vconstraint,
        mapping.msgs_names,
        messages,
        mapping.marginals_names,
        marginals,
        mapping.meta,
        mapping.factornode
    )

    return Message(message, is_message_clamped, is_message_initial)
end

Base.map(::Type{T}, mapping::M) where { T, M <: MessageMapping } = Rocket.MapOperator{T, M}(mapping)


text:
import Base.isequal
import Base.hash
import Base.show
import Base.isless

"""
EMURxn has fields:
    flux::String, product_emu::ProductEMU, reactand_emus::Vector{ReactandEMU}

and represents an EMU Reaction
"""
struct EMURxn
    flux::String
    product_emu::ProductEMU
    reactand_emus::Vector{ReactandEMU}
    #transition::Dict{RAS,PAS}
end

"""
SpeciesEMURxn has fields:
    flux::String, product_emu::SpeciesEMU, reactand_emus::Vector{SpeciesEMU}

and represents an EMU Reaction on the species / SAS level
"""
struct SpeciesEMURxn
    flux::String
    product_emu::SpeciesEMU
    reactand_emus::Vector{SpeciesEMU}
    #transition::Dict{RAS,PAS}
end

# order: 1. product emu, 2. flux , this should be ok
function Base.isless(a::SpeciesEMURxn,b::SpeciesEMURxn)
    if(a.product_emu<b.product_emu) ; return true; end
    if(a.product_emu==b.product_emu)
        return a.flux < b.flux
    else
        return false
    end
end
function Base.isequal(a::SpeciesEMURxn,b::SpeciesEMURxn)
    if(a.flux < b.flux )
        return true
    end
    if(f.flux == b.flux)
        if( isequal(a.product_emu , b.product_emu) && a.flux ==b.flux )
            if(length(a.reactand_emus)!=lenght(b.reactand_emus))
                return false
            else
                rea = sort(a.reactand_emus)
                reb = sort(b.reactand_emus)
                equal = all( isequal.(rea,reb) )
                return equal
            end
        else
            return false
        end
    else
        return false
    end
    return false
end
function Base.hash(x::SpeciesEMURxn)
    hash( string("semurxn$(x.flux):$(hash(x.Set(product_emu))):$(hash(Set(x.educt_emus)))" )  )
end

function SpeciesEMURxn(r::EMURxn)
    SpeciesEMURxn( r.flux , pemu2semu( r.product_emu ) , map( x1 -> remu2semu(x1) , r.reactand_emus ) )
end
function Base.show(io::IO,r::SpeciesEMURxn)
    print(io,"$(r.flux): ")
    for si in r.reactand_emus
        print(io, "$(si) " )
    end
    print(io," --> ")
    print(io,r.product_emu)
end

# compares the flux and the product emu, this should uniquely identify the emurxn
function Base.isequal(r1::EMURxn,r2::EMURxn)
    return r1.flux==r2.flux && r1.product_emu==r2.product_emu
end
function Base.hash(r::EMURxn)
    return hash( string(r.flux, "$(hash(r.product_emu))" ) )
end

"""
compute_emu_decomp(n::Network,targets::Vector{Set{SAS}})
returns Vector{EMURxn}

computes EMU decomposition for the given network and the given target SpeciesEMUs
"""
function compute_emu_decomp(n::Network,targets::Vector{Set{SAS}})

    debuglevel=1

    # prepare:
    # compute mappings SAS <-> RAS , SAS <-> PAS
    all_sas::Set{SAS} = sas(n)
    all_ras::Set{RAS} = ras(n)
    all_pas::Set{PAS} = pas(n)
    d_sas2pas::Dict{SAS,Set{PAS}} = Dict{SAS,Set{PAS}}()
    d_pas2sas::Dict{PAS,SAS}      = Dict{PAS,SAS}()
    d_sas2ras::Dict{SAS,Set{RAS}} = Dict{SAS,Set{RAS}}()
    d_ras2sas::Dict{RAS,SAS}      = Dict{RAS,SAS}()
    for si in all_sas
        m_ras::Set{RAS} = filter( x1 -> x1.sas==si , all_ras )
        d_sas2ras[si] = Set(m_ras)
        for ri in m_ras
            d_ras2sas[ri] = si
        end
        m_pas::Set{PAS} = filter( x1 -> x1.sas==si , all_pas )
        d_sas2pas[si] = Set(m_pas)
        for pi in m_pas
            d_pas2sas[pi] = si
        end
    end

    if(debuglevel>=2)
        println("\n------------------------------------------------------\n")
        println(d_sas2pas)
        println("\n------------------------------------------------------\n")
        println(d_sas2ras)
        println("\n")
    end

    ## Start with decomposition:

    processed_emus::Set{SpeciesEMU} = Set{SpeciesEMU}()
    active_emus::Set{SpeciesEMU} = Set{SpeciesEMU}()

    ## check if we have target semus which span more than one metabolite, those will be split..
    for ti in targets
        if(!isinsinglespecies(ti))
            printf("Encountered multi-species target SEMUs, those will be split..")
            tis = splitSEMUIntoSpecies(target)
            for tisi in values(tis)
                push!( active_emus , tisi )
            end
        else
            push!( active_emus , ti )
        end
    end

    emu_rxns::Vector{EMURxn} = Vector{EMURxn}();

    while( ! isempty(active_emus) )
        active_emus = setdiff(active_emus,processed_emus)
        if(isempty(active_emus))
            break # we' 're done
        end

        # get next emu
        emu = pop!(active_emus)
        push!(processed_emus,emu)

        if(debuglevel>=2); println(string("Process EMU ",emu)); end

        # 1. construct EMU Reactions which gives rise to this SEMU
        # 1.1 find all PAS which contain this Species EMU
        pemus::Set{PAS} = findPASsForSEMU( d_sas2pas , emu )
        # 1.2 split into "products", product is (flux,species,ppos)
        sorted_pemus::Dict{Tuple{String,String,Int64},Set{PAS}} = splitPASsIntoProducts(pemus)

        # 2 loop over sorted_pemus and construct emu rxn for each

        for pemu_p::Tuple{String,String,Int64} in keys(sorted_pemus)
            if( isinflux( n.fluxes_db[pemu_p[1]] ) )
                if(debuglevel>=2); println("Hit Influx $(pemu_p[1])"); end
                continue
            end

            pemu::Set{PAS} = sorted_pemus[pemu_p]
            flux_i = n.fluxes_db[pemu_p[1]]
            rtransition = map(reverse, flux_i.transition)

            # 2.1 find all ras:
            # sort pemu!..
            vpemu::Vector{PAS} = sort(collect(pemu))
            vras::Vector{RAS} = map( x1 -> rtransition[x1] , vpemu )
            # 1.3 Split REMUs into reactands, "reactand" is (flux,species,rpos)
            sorted_ras::Dict{Tuple{String,String,Int64},Set{RAS}} = splitRASsIntoReactands(Set(vras))

            # create EMU Rxn..
            if(debuglevel>=2)
                println("Make EMU Rxn:")
                println("Target: \n $(vpemu)")
                println("REMUs Sizes: $(map( x -> length(sorted_ras[x]),keys(sorted_ras))) ")
                println("REMUs:")
                for ki in keys(sorted_ras)
                    println(sorted_ras[ki])
                end
                # \n  $(sorted_ras)\n")
            end

            # create rxn
            push!(emu_rxns, EMURxn( pemu_p[1] , pemu , collect(values(sorted_ras)) ) )

            # add mew semus to be processed:
            for ki in keys(sorted_ras)
                push!( active_emus , remu2semu(sorted_ras[ki]) )
            end
        end
    end
    return emu_rxns
end

#return Dict{String,Set{SAS}}
function splitSEMUIntoSpecies(sass::Set{SAS})
    all_species = map( x::SAS -> (x.species) , sass )
    species = Set(all_species)
    sortedspecies = Dict{String,Set{SAS}}()
    for si in species
        sortedspecies[si] = Set( filter( x1 -> x1.species==si , sass) )
    end
    return sortedspecies
end

function splitPASsIntoFluxes(sp::Set{PAS})
    pfluxes = Set(map( x::PAS -> x.flux , pemus ))
    sortedpas = Dict{String,Set{PAS}}()
    for pi in pfluxes
        sortedpas[si] = Set( filter( x1 -> x1.flux==pi , sp) )
    end
    return sortedpas
end

# product is ( "Flux" , "Species" , "ppos" )
function splitPASsIntoProducts(sp::Set{PAS})
    products::Set{Tuple{String,String,Int64}} =  Set{Tuple{String,String,Int64}}( map( x -> (x.flux,x.sas.species,x.ppos) , sp ) )
    sorted_pas = Dict{Tuple{String,String,Int64},Set{PAS}}()
    for pi in products
        sorted_pas[pi] = Set( filter( x1 -> ( (x1.flux==pi[1]) && (x1.sas.species==pi[2]) && (x1.ppos==pi[3]) ), sp) )
    end
    return sorted_pas
end

function isinsingleflux(sp::Set{RAS})
    if(isempty(sp)); return true; end
    sp1 = collect(sp)[1].flux
    return all( map( x -> x.flux == sp1 , sp ) )
end

# ! Function assumes that all are in same flux !
# reactand is ( "Flux" , "Species" , "rpos" )
function splitRASsIntoReactands(sp::Set{RAS})
    if(!isinsingleflux(sp)); error("::findProductEMUsForSEMU : predondition violated"); end
    reactands::Set{Tuple{String,String,Int64}} =  Set{Tuple{String,String,Int64}}( map( x -> (x.flux,x.sas.species,x.rpos) , sp ) )
    sorted_ras::Dict{Tuple{String,String,Int64},Set{RAS}} = Dict{Tuple{String,String,Int64},Set{RAS}}()
    for ri in reactands
        sorted_ras[ri] = Set( filter( x1 -> ( (x1.flux==ri[1]) && (x1.sas.species==ri[2]) && (x1.rpos==ri[3]) ), sp) )
    end
    return sorted_ras
end

# assumption: emu is inside one species!
function findPASsForSEMU( d_sas2pas::Dict{SAS,Set{PAS}} , emu::SpeciesEMU )
    # sanity check:
    if(!isinsinglespecies(emu)); error("::findProductEMUsForSEMU : predondition violated"); end

    # all PAS :
    all_pas::Set{PAS} = reduce( (a,b) -> union(a,b) , map( x -> d_sas2pas[x] , emu ) )

    return all_pas
end


text:
include("./Measurements.jl")

abstract type AbstractCovarianceFunction{Tv} end

# automatically implements a mutating batched version of a given covariance function
function (cov::AbstractCovarianceFunction{Tv})(out::AbstractMatrix{Tv}, x_vec::AbstractVector{<:AbstractMeasurement}, y_vec::AbstractVector{<:AbstractMeasurement}) where Tv
    for cartesian in CartesianIndices(out)
        out[cartesian] = cov(x_vec[cartesian[1]], y_vec[cartesian[2]])
    end
end

# automatically implements a mutating batched version of a given covariance function, using symmetry
function (cov::AbstractCovarianceFunction{Tv})(out::AbstractMatrix{Tv}, x_vec::AbstractVector{<:AbstractMeasurement}) where Tv
    for cartesian in CartesianIndices(out)
        if cartesian[1] >= cartesian[2]
            out[cartesian] = cov(x_vec[cartesian[1]], x_vec[cartesian[2]])
        else
            out[cartesian] = out[cartesian[2], cartesian[1]]
        end
    end
end

struct MaternCovariance1_2{Tv}<:AbstractCovarianceFunction{Tv}
    length_scale::Tv
end
function (cov::MaternCovariance1_2)(x::PointMeasurement, y::PointMeasurement)
    dist = norm(x.coordinate - y.coordinate);
    sigma = cov.length_scale;
    return exp(-dist/sigma)
end

struct MaternCovariance3_2{Tv}<:AbstractCovarianceFunction{Tv}
    length_scale::Tv
end
function (cov::MaternCovariance3_2)(x::PointMeasurement, y::PointMeasurement)
    dist = norm(x.coordinate - y.coordinate);
    sigma = cov.length_scale;
    return (1+sqrt(3)*dist/sigma)*exp(-sqrt(3)*dist/sigma)
end

struct MaternCovariance5_2{Tv}<:AbstractCovarianceFunction{Tv}
    length_scale::Tv
end

function (cov::MaternCovariance5_2)(x::PointMeasurement, y::PointMeasurement)
    # d = length(x.coordinate)
    dist = norm(x.coordinate - y.coordinate);
    sigma = cov.length_scale;
    F(t,a) = (1 + sqrt(5)*t/a + 5*t^2/(3*a^2))*exp(-sqrt(5)*t/a)
    return F(dist,sigma)
end

function (cov::MaternCovariance5_2)(x::ΔδPointMeasurement, y::ΔδPointMeasurement)
    d = length(x.coordinate);
    w1_x = x.weight_Δ;
    w2_x = x.weight_δ;
    w1_y = y.weight_Δ;
    w2_y = y.weight_δ;

    F(t,a) = (1 + sqrt(5)*t/a + 5*t^2/(3*a^2))*exp(-sqrt(5)*t/a)
    D2F(t,a) = -5*(d*a^2+sqrt(5)*d*a*t-5*t^2)/(3*a^4) * exp(-sqrt(5)*t/a);
    D4F(t,a) = 25*(d*(d+2)*a^2-(3+2*d)*sqrt(5)*a*t+5*t^2)/(3*a^6) * exp(-sqrt(5)*t/a);
    dist = norm(x.coordinate - y.coordinate);
    sigma = cov.length_scale;
    return w1_x*w1_y*D4F(dist,sigma) + (w2_x*w1_y+w1_x*w2_y)*D2F(dist,sigma) + w2_x*w2_y*F(dist,sigma)
end

function (cov::MaternCovariance5_2)(x::Δ∇δPointMeasurement, y::Δ∇δPointMeasurement)
    d = length(x.coordinate);
    w1_x = x.weight_Δ;
    w2_x = x.weight_δ;
    wg_x = x.weight_∇
    w1_y = y.weight_Δ;
    w2_y = y.weight_δ;
    wg_y = y.weight_∇

    F(t,a) = (1 + sqrt(5)*t/a + 5*t^2/(3*a^2))*exp(-sqrt(5)*t/a)
    D2F(t,a) = -5*(d*a^2+sqrt(5)*d*a*t-5*t^2)/(3*a^4) * exp(-sqrt(5)*t/a)
    D4F(t,a) = 25*(d*(d+2)*a^2-(3+2*d)*sqrt(5)*a*t+5*t^2)/(3*a^6) * exp(-sqrt(5)*t/a)
    DF(t,a) = -5*(a+sqrt(5)*t)*exp(-sqrt(5)*t/a)/(3*a^3)
    D3F(t,a) = 25*exp(-sqrt(5)*t/a)*(a*(2+d)-sqrt(5)*t)/(3*a^5)
    DDF(t,a) = 25*exp(-sqrt(5)*t/a)/(3*a^4)
    vec = x.coordinate - y.coordinate
    dist = norm(vec);
    sigma = cov.length_scale;
    return w1_x*w1_y*D4F(dist,sigma) + (w2_x*w1_y+w1_x*w2_y)*D2F(dist,sigma) + w2_x*w2_y*F(dist,sigma) - w1_x*D3F(dist,sigma)*sum(vec.*wg_y) + w1_y*D3F(dist,sigma)*sum(vec.*wg_x) - w2_x*DF(dist,sigma)*sum(vec.*wg_y) + w2_y*DF(dist,sigma)*sum(vec.*wg_x) + (sum(-wg_x.*wg_y)*DF(dist,sigma)+sum(wg_x.*vec)*sum(-wg_y.*vec)*DDF(dist,sigma))
end

struct MaternCovariance7_2{Tv}<:AbstractCovarianceFunction{Tv}
    length_scale::Tv
end

function (cov::MaternCovariance7_2)(x::PointMeasurement, y::PointMeasurement)
    dist = norm(x.coordinate - y.coordinate);
    sigma = cov.length_scale;
    F(t,a) = (15*a^3+15*sqrt(7)*a^2*t+42*a*t^2+7*sqrt(7)*t^3)/(15*a^3)*exp(-sqrt(7)*t/a);
    return F(dist,sigma)
end

function (cov::MaternCovariance7_2)(x::ΔδPointMeasurement, y::ΔδPointMeasurement)
    d = length(x.coordinate);
    w1_x = x.weight_Δ;
    w2_x = x.weight_δ;
    w1_y = y.weight_Δ;
    w2_y = y.weight_δ;
    F(t,a) = (15*a^3+15*sqrt(7)*a^2*t+42*a*t^2+7*sqrt(7)*t^3)/(15*a^3)*exp(-sqrt(7)*t/a);
    D2F(t,a) = -7*(3*d*a^3+3*sqrt(7)*a^2*d*t+7*a*(d-1)*t^2-7*sqrt(7)*t^3)/(15*a^5)*exp(-sqrt(7)*t/a);
    D4F(t,a) = 49*(d*(d+2)*a^3+d*(d+2)*sqrt(7)*a^2*t-14*a*(2+d)*t^2+7*sqrt(7)*t^3)/(15*a^7)*exp(-sqrt(7)*t/a);
    dist = norm(x.coordinate - y.coordinate);
    sigma = cov.length_scale;
    return w1_x*w1_y*D4F(dist,sigma) + (w2_x*w1_y+w1_x*w2_y)*D2F(dist,sigma) + w2_x*w2_y*F(dist,sigma)
end

function (cov::MaternCovariance7_2)(x::Δ∇δPointMeasurement, y::Δ∇δPointMeasurement)
    d = length(x.coordinate);
    w1_x = x.weight_Δ;
    w2_x = x.weight_δ;
    wg_x = x.weight_∇
    w1_y = y.weight_Δ;
    w2_y = y.weight_δ;
    wg_y = y.weight_∇

    F(t,a) = (15*a^3+15*sqrt(7)*a^2*t+42*a*t^2+7*sqrt(7)*t^3)/(15*a^3)*exp(-sqrt(7)*t/a);
    D2F(t,a) = -7*(3*d*a^3+3*sqrt(7)*a^2*d*t+7*a*(d-1)*t^2-7*sqrt(7)*t^3)/(15*a^5)*exp(-sqrt(7)*t/a);
    D4F(t,a) = 49*(d*(d+2)*a^3+d*(d+2)*sqrt(7)*a^2*t-14*a*(2+d)*t^2+7*sqrt(7)*t^3)/(15*a^7)*exp(-sqrt(7)*t/a);
    DF(t,a) = -7*(3*a^2+3*sqrt(7)*a*t+7*t^2)*exp(-sqrt(7)*t/a)/(15*a^4)
    D3F(t,a) = 49*exp(-sqrt(7)*t/a)*(a^2*(2+d)+sqrt(7)*a*(2+d)*t-7*t^2)/(15*a^6)
    DDF(t,a) = 49*exp(-sqrt(7)*t/a)*(a+sqrt(7)*t)/(15*a^5)
    vec = x.coordinate - y.coordinate
    dist = norm(vec);
    sigma = cov.length_scale;
    return w1_x*w1_y*D4F(dist,sigma) + (w2_x*w1_y+w1_x*w2_y)*D2F(dist,sigma) + w2_x*w2_y*F(dist,sigma) - w1_x*D3F(dist,sigma)*sum(vec.*wg_y) + w1_y*D3F(dist,sigma)*sum(vec.*wg_x) - w2_x*DF(dist,sigma)*sum(vec.*wg_y) + w2_y*DF(dist,sigma)*sum(vec.*wg_x) + (sum(-wg_x.*wg_y)*DF(dist,sigma)+sum(wg_x.*vec)*sum(-wg_y.*vec)*DDF(dist,sigma))
end

struct MaternCovariance9_2{Tv}<:AbstractCovarianceFunction{Tv}
    length_scale::Tv
end

function (cov::MaternCovariance9_2)(x::PointMeasurement, y::PointMeasurement)
    dist = norm(x.coordinate - y.coordinate);
    sigma = cov.length_scale;
    F(t,a) = (35*a^4+105*a^3*t+135*a^2*t^2+90*a*t^3+27*t^4)/(35*a^4)*exp(-3*t/a);
    return F(dist,sigma)
end

function (cov::MaternCovariance9_2)(x::ΔδPointMeasurement, y::ΔδPointMeasurement)
    d = length(x.coordinate);
    w1_x = x.weight_Δ;
    w2_x = x.weight_δ;
    w1_y = y.weight_Δ;
    w2_y = y.weight_δ;
    F(t,a) = (35*a^4+105*a^3*t+135*a^2*t^2+90*a*t^3+27*t^4)/(35*a^4)*exp(-3*t/a);
    D2F(t,a) = -9*(5*d*a^4+15*d*a^3*t+9*a^2*(2*d-1)*t^2+9*a*(d-3)*t^3-27*t^4)/(35*a^6)*exp(-3*t/a);
    D4F(t,a) = 81*(d*(d+2)*a^4+3*a^3*d*(d+2)*t+3*a^2*(d^2-4)-18*a*(d+2)*t^3+27*t^4)/(35*a^8)*exp(-3*t/a);
    dist = norm(x.coordinate - y.coordinate);
    sigma = cov.length_scale;
    return w1_x*w1_y*D4F(dist,sigma) + (w2_x*w1_y+w1_x*w2_y)*D2F(dist,sigma) + w2_x*w2_y*F(dist,sigma)
end

function (cov::MaternCovariance9_2)(x::Δ∇δPointMeasurement, y::Δ∇δPointMeasurement)
    d = length(x.coordinate);
    w1_x = x.weight_Δ;
    w2_x = x.weight_δ;
    wg_x = x.weight_∇
    w1_y = y.weight_Δ;
    w2_y = y.weight_δ;
    wg_y = y.weight_∇

    F(t,a) = (35*a^4+105*a^3*t+135*a^2*t^2+90*a*t^3+27*t^4)/(35*a^4)*exp(-3*t/a);
    D2F(t,a) = -9*(5*d*a^4+15*d*a^3*t+9*a^2*(2*d-1)*t^2+9*a*(d-3)*t^3-27*t^4)/(35*a^6)*exp(-3*t/a);
    D4F(t,a) = 81*(d*(d+2)*a^4+3*a^3*d*(d+2)*t+3*a^2*(d^2-4)*(t^2)-18*a*(d+2)*t^3+27*t^4)/(35*a^8)*exp(-3*t/a);
    DF(t,a) = -9*(5*a^3+15*a^2*t+18*a*t^2+9*t^3)*exp(-3*t/a)/(35*a^5)
    D3F(t,a) = 81*exp(-3*t/a)*(a^3*(2+d)+3*a^2*(2+d)*t+3*a*(1+d)*t^2-9*t^3)/(35*a^7)
    DDF(t,a) = 81*exp(-3*t/a)*(a^2+3*a*t+3*t^2)/(35*a^6)
    vec = x.coordinate - y.coordinate
    dist = norm(vec);
    sigma = cov.length_scale;
    return w1_x*w1_y*D4F(dist,sigma) + (w2_x*w1_y+w1_x*w2_y)*D2F(dist,sigma) + w2_x*w2_y*F(dist,sigma) - w1_x*D3F(dist,sigma)*sum(vec.*wg_y) + w1_y*D3F(dist,sigma)*sum(vec.*wg_x) - w2_x*DF(dist,sigma)*sum(vec.*wg_y) + w2_y*DF(dist,sigma)*sum(vec.*wg_x) + (sum(-wg_x.*wg_y)*DF(dist,sigma)+sum(wg_x.*vec)*sum(-wg_y.*vec)*DDF(dist,sigma))
end

struct MaternCovariance11_2{Tv}<:AbstractCovarianceFunction{Tv}
    length_scale::Tv
end

# Matern covariance function
function (cov::MaternCovariance11_2)(x::PointMeasurement, y::PointMeasurement)
    dist = norm(x.coordinate - y.coordinate);
    sigma = cov.length_scale;
    F(t,a) = (945*a^5+945*sqrt(11)*a^4*t+4620*a^3*t^2+1155*sqrt(11)*a^2*t^3+1815*a*t^4+121*sqrt(11)*t^5)/(945*a^5)*exp(-sqrt(11)*t/a);
    return F(dist,sigma)
end

function (cov::MaternCovariance11_2)(x::ΔδPointMeasurement, y::ΔδPointMeasurement)
    d = length(x.coordinate);
    w1_x = x.weight_Δ;
    w2_x = x.weight_δ;
    w1_y = y.weight_Δ;
    w2_y = y.weight_δ;
    F(t,a) = (945*a^5+945*sqrt(11)*a^4*t+4620*a^3*t^2+1155*sqrt(11)*a^2*t^3+1815*a*t^4+121*sqrt(11)*t^5)/(945*a^5)*exp(-sqrt(11)*t/a);
    D2F(t,a) = -11*(105*d*a^5+105*d*sqrt(11)*a^4*t+165*(3*d-1)*a^3*t^2+55*sqrt(11)*(2*d-3)*a^2*t^3+121*(d-6)*a*t^4-121*sqrt(11)*t^5)/(945*a^7)*exp(-sqrt(11)*t/a);
    D4F(t,a) = 121*(15*d*(2+d)*a^5+15*d*(2+d)*sqrt(11)*a^4*t+66*(d^2+d-2)*a^3*t^2+11*sqrt(11)*a^2*(d^2-4*d-12)*t^3-121*(3+2*d)*a*t^4+121*sqrt(11)*t^5)/(945*a^9)*exp(-sqrt(11)*t/a);
    dist = norm(x.coordinate - y.coordinate);
    sigma = cov.length_scale;
    return w1_x*w1_y*D4F(dist,sigma) + (w2_x*w1_y+w1_x*w2_y)*D2F(dist,sigma) + w2_x*w2_y*F(dist,sigma)
end

function (cov::MaternCovariance11_2)(x::Δ∇δPointMeasurement, y::Δ∇δPointMeasurement)
    d = length(x.coordinate);
    w1_x = x.weight_Δ;
    w2_x = x.weight_δ;
    wg_x = x.weight_∇
    w1_y = y.weight_Δ;
    w2_y = y.weight_δ;
    wg_y = y.weight_∇

    F(t,a) = (945*a^5+945*sqrt(11)*a^4*t+4620*a^3*t^2+1155*sqrt(11)*a^2*t^3+1815*a*t^4+121*sqrt(11)*t^5)/(945*a^5)*exp(-sqrt(11)*t/a);
    D2F(t,a) = -11*(105*d*a^5+105*d*sqrt(11)*a^4*t+165*(3*d-1)*a^3*t^2+55*sqrt(11)*(2*d-3)*a^2*t^3+121*(d-6)*a*t^4-121*sqrt(11)*t^5)/(945*a^7)*exp(-sqrt(11)*t/a);
    D4F(t,a) = 121*(15*d*(2+d)*a^5+15*d*(2+d)*sqrt(11)*a^4*t+66*(d^2+d-2)*a^3*t^2+11*sqrt(11)*a^2*(d^2-4*d-12)*t^3-121*(3+2*d)*a*t^4+121*sqrt(11)*t^5)/(945*a^9)*exp(-sqrt(11)*t/a);
    DF(t,a) = -11*(105*a^4+105*sqrt(11)*a^3*t+495*a^2*t^2+110*sqrt(11)*a*t^3+121*t^4)*exp(-sqrt(11)*t/a)/(945*a^6)
    D3F(t,a) = 121*exp(-sqrt(11)*t/a)*(15*a^4*(2+d)+15*sqrt(11)*a^3*(2+d)*t+33*a^2*(3+2*d)*t^2+11*sqrt*(11)*a*(d-1)*t^3-121*t^4)/(945*a^8)
    DDF(t,a) = 121*exp(-sqrt(11)*t/a)*(15*a^3+15*sqrt(11)*a^2*t+66*a*t^2+11*sqrt(11)*t^3)/(945*a^7)
    vec = x.coordinate - y.coordinate
    dist = norm(vec);
    sigma = cov.length_scale;
    return w1_x*w1_y*D4F(dist,sigma) + (w2_x*w1_y+w1_x*w2_y)*D2F(dist,sigma) + w2_x*w2_y*F(dist,sigma) - w1_x*D3F(dist,sigma)*sum(vec.*wg_y) + w1_y*D3F(dist,sigma)*sum(vec.*wg_x) - w2_x*DF(dist,sigma)*sum(vec.*wg_y) + w2_y*DF(dist,sigma)*sum(vec.*wg_x) + (sum(-wg_x.*wg_y)*DF(dist,sigma)+sum(wg_x.*vec)*sum(-wg_y.*vec)*DDF(dist,sigma))
end

# The Gaussian covariance function
struct GaussianCovariance{Tv}<:AbstractCovarianceFunction{Tv}
    length_scale::Tv
end

function (cov::GaussianCovariance)(x::PointMeasurement, y::PointMeasurement)
    dist = norm(x.coordinate - y.coordinate);
    sigma = cov.length_scale;
    return exp(-dist^2/(2*sigma^2))
end

function (cov::GaussianCovariance)(x::ΔδPointMeasurement, y::ΔδPointMeasurement)
    d = length(x.coordinate);
    w1_x = x.weight_Δ;
    w2_x = x.weight_δ;
    w1_y = y.weight_Δ;
    w2_y = y.weight_δ;
    F(t,a) = exp(-t^2/(2*a^2))
    D2F(t,a) = (t^2 - a^2*d)/(a^4)*exp(-t^2/(2*a^2));
    D4F(t,a) = (a^4*d*(2+d)-2*a^2*(2+d)*t^2+t^4)*exp(-t^2/(2*a^2))/a^8
    dist = norm(x.coordinate - y.coordinate);
    sigma = cov.length_scale;
    return w1_x*w1_y*D4F(dist,sigma) + (w2_x*w1_y+w1_x*w2_y)*D2F(dist,sigma) + w2_x*w2_y*F(dist,sigma)
end

function (cov::GaussianCovariance)(x::Δ∇δPointMeasurement, y::Δ∇δPointMeasurement)
    d = length(x.coordinate);
    w1_x = x.weight_Δ;
    w2_x = x.weight_δ;
    wg_x = x.weight_∇
    w1_y = y.weight_Δ;
    w2_y = y.weight_δ;
    wg_y = y.weight_∇

    F(t,a) = exp(-t^2/(2*a^2))
    D2F(t,a) = (t^2 - a^2*d)/(a^4)*exp(-t^2/(2*a^2));
    D4F(t,a) = (a^4*d*(2+d)-2*a^2*(2+d)*t^2+t^4)*exp(-t^2/(2*a^2))/a^8
    DF(t,a) = -exp(-t^2/(2*a^2))/a^2
    D3F(t,a) = exp(-t^2/(2*a^2))*(a^2*(2+d)*t-t^3)/a^6
    DDF(t,a) = exp(-t^2/(2*a^2))/a^4
    vec = x.coordinate - y.coordinate
    dist = norm(vec);
    sigma = cov.length_scale;
    return w1_x*w1_y*D4F(dist,sigma) + (w2_x*w1_y+w1_x*w2_y)*D2F(dist,sigma) + w2_x*w2_y*F(dist,sigma) - w1_x*D3F(dist,sigma)*sum(vec.*wg_y) + w1_y*D3F(dist,sigma)*sum(vec.*wg_x) - w2_x*DF(dist,sigma)*sum(vec.*wg_y) + w2_y*DF(dist,sigma)*sum(vec.*wg_x) + (sum(-wg_x.*wg_y)*DF(dist,sigma)+sum(wg_x.*vec)*sum(-wg_y.*vec)*DDF(dist,sigma))
end

function (cov::AbstractCovarianceFunction)(x::ΔδPointMeasurement, y::PointMeasurement)
    return cov(x, ΔδPointMeasurement(y))
end
function (cov::AbstractCovarianceFunction)(x::PointMeasurement, y::ΔδPointMeasurement)
    return (cov::AbstractCovarianceFunction)(y,x)
end

function (cov::AbstractCovarianceFunction)(x::Δ∇δPointMeasurement, y::PointMeasurement)
    return cov(x, Δ∇δPointMeasurement(y))
end
function (cov::AbstractCovarianceFunction)(x::PointMeasurement, y::Δ∇δPointMeasurement)
    return (cov::AbstractCovarianceFunction)(y,x)
end


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

include("geoops/filter.jl")
include("geoops/groupby.jl")
include("geoops/unique.jl")
include("geoops/integrate.jl")


text:
module Flutter

function dist(x)
    [sqrt(x[i,1]*x[i,1]+x[i,2]*x[i,2]+1e4*x[i,3]*x[i,3]) for i=1:size(x,1)]
end

function read(path)
    readdlm(path,'\t';header=true)
end

function readdist(path)
    x,h=read(path)
    [x[:,1] dist(x[:,2:4]) dist(x[:,5:7]) dist(x[:,8:10])],h
end

function diff(x)
    [x[2:end,1] diff(x[:,2:end])]
end

end

max_stars_repo_name: jlep/flutter

max_issues_repo_name: jlep/flutter

max_forks_repo_name: jlep/flutter


text:
export DegreeVertexCover

struct DegreeVertexCover end

"""
    vertex_cover(g, DegreeVertexCover())

Obtain a vertex cover using a greedy heuristic.

### Implementation Notes
An edge is said to be covered if it has at least one end-point in the vertex cover.
Initialise the vertex cover to an empty set and iteratively choose the vertex with the most uncovered
edges.

### Performance
Runtime: O((|V|+|E|)*log(|V|))
Memory: O(|V|)

# Examples
```jldoctest
julia> using LightGraphs

julia> vertex_cover(PathGraph(3), DegreeVertexCover())
1-element Array{Int64,1}:
 2

julia> vertex_cover(CycleGraph(3), DegreeVertexCover())
2-element Array{Int64,1}:
 1
 3
```
"""
function vertex_cover(
    g::AbstractGraph{T},
    alg::DegreeVertexCover
    ) where T <: Integer

    nvg = nv(g)
    in_cover = falses(nvg)
    length_cover = 0
    degree_queue = PriorityQueue(Base.Order.Reverse, enumerate(degree(g)))

    while !isempty(degree_queue) && peek(degree_queue)[2] > 0
        v = dequeue!(degree_queue)
        in_cover[v] = true
        length_cover += 1

        @inbounds @simd for u in neighbors(g, v)
            if !in_cover[u]
                degree_queue[u] -= 1
            end
        end
    end
    return LightGraphs.findall!(in_cover, Vector{T}(undef, length_cover))
end


text:
using Sockets

server = Sockets.listen(8080)

while true
    sock = Sockets.accept(server)
    @async begin
        data = readline(sock)
        print("Got request:\n", data, "\n")
        cmd = split(data, " ")[2][2:end]
        println(sock, "\nHTTP/1.1 200 OK\nContent-Type: text/html\n")
        println(sock, string("", cmd, "=", eval(Meta.parse(cmd)), ""))
        close(sock)
    end
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using BangBang
using LinearAlgebra

function rnn!!(n, J, x)
    dest = x
    y = similar(x)
    for _ in 1:n
        @! y = mul!(y, J, x)
        @! y .= tanh.(y)
        x, y = y, x
    end
    @! dest .= x
    return dest
end

using Zygote
d = 10
J = randn(d, d)
x0 = randn(d)
y_target = randn(d)
g, = Zygote.gradient(J -> sum((rnn!!(20, J, x0) .- y_target) .^ 2), J)


text:
using TKPappKitchenSink
using Base.Test
using StaticArrays

@testset "similar basis vector" begin
    @test similar_basis_vector(randn(3), 1) == [1.0, 0.0, 0.0]
    @test similar_basis_vector(randn(3), 2) == [0.0, 1.0, 0.0]
    @test similar_basis_vector(SVector(1.0, 2.0), 1) ≡ SVector(1.0, 0.0)
    @test similar_basis_vector(SVector(1.0, 2.0), 2) ≡ SVector(0.0, 1.0)
end

@testset "function slice" begin
    f(x) = sum(abs2, x)
    x = randn(3)
    @test function_slice(f, x, 1)(0.3) == f(x + 0.3*[1.0, 0.0, 0.0])
end


text:
# loading the Julia packages  needed.
using DifferentialEquations
using LinearAlgebra, DiffEqSensitivity, Optim
using Flux: flatten, params, relu
using DiffEqFlux, Flux
using Plots
using Flux: train!
using GalacticOptim
using Optim
using Turing, Distributions
using MCMCChains, StatsPlots
using CSV, DataFrames
using SymbolicRegression
using Random
Random.seed!(14);
source_data = DataFrame(CSV.File("./DeepLearningEffectiveReproductionNumber/Source_Data/Provincial_Daily_Totals.csv"))
data_on = source_data[source_data.Province.=="ONTARIO", :]
n = 30
m = 149
data_acc = data_on.TotalCases[(n+1):n+m+1]
data_daily = data_on.TotalCases[(n+1):n+m+1] - data_on.TotalCases[n:n+m]
display(plot(data_daily, label = "Daily Confirmed Cases", lw = 2))
display(plot(data_acc, label = "Accumulated Confirmed Cases", lw = 2))
data_daily[1]
println(length(data_acc))

# Model generation
ann = FastChain(FastDense(1, 32, tanh), FastDense(32, 1))
p_0 = initial_params(ann)
function SIR_nn(du, u, p, t)
    I, H = u
    du[1] = 0.1 * min(5, abs(ann(t, p)[1])) * I - 0.1 * I
    du[2] = 0.1 * min(5, abs(ann(t, p)[1])) * I
end
u_0 = Float32[1, data_acc[1]]
tspan_data = (0.0f0, 149.0f0)
prob_nn = ODEProblem(SIR_nn, u_0, tspan_data, p_0)
function train(θ)
    solvedata = Array(concrete_solve(prob_nn, Vern7(), u_0, θ, saveat = 1,
        abstol = 1e-6, reltol = 1e-6,
        sensealg = InterpolatingAdjoint(autojacvec = ReverseDiffVJP())))
end
function loss(θ)
    pred = train(θ)
    sum(abs2, (log.(data_acc) .- log.(pred[2, 1:150]))), pred # + 1e-5*sum(sum.(abs, params(ann)))
end

const losses = []
callback(θ, l, pred) = begin
    push!(losses, l)
    if length(losses) % 50 == 0
        println(losses[end])
    end
    false
end

@time res1_node = DiffEqFlux.sciml_train(loss, p_0, ADAM(0.06), cb = callback, maxiters = 1000)
@time res2_node = DiffEqFlux.sciml_train(loss, res1_node.minimizer, BFGS(initial_stepnorm = 0.01), cb = callback, maxiters = 300)
ann_param = res2_node.minimizer

using BSON: @save
@save "./DeepLearningEffectiveReproductionNumber/Saving_Data/ann_nn_ir.bason" ann
@save "./DeepLearningEffectiveReproductionNumber/Saving_Data/ann_para_irbfgs100.bason" ann_param
@time res3_node = DiffEqFlux.sciml_train(loss, res2_node.minimizer, BFGS(initial_stepnorm = 0.01), cb = callback, maxiters = 200)
ann_param = res3_node.minimizer
@save "./DeepLearningEffectiveReproductionNumber/Saving_Data/ann_para_irbfgs500.bason" ann_param

p_min = ann_param
tspan_predict = (0.0, 149)
scatter(data_acc, label = "Real accumulated cases")
prob_prediction = ODEProblem(SIR_nn, u_0, tspan_data, p_min)
data_prediction = Array(solve(prob_prediction, Tsit5(), saveat = 1))
plot!(data_prediction[2, :], label = "Fit accumulated cases")

xlabel!("Days after Feb 25")

tspan = collect(0:1:149)'
ann_value = abs.(ann(tspan, p_min))'
plot(ann_value, label = "Effective reproduction number")
xlabel!("Days after Feb 25")


text:
using Printf

struct Sphere
    sphere_id::Int64
    cluster_id::Int64
    coordinates::Array{Float64, 1}
    line::String
end

function GetClusterHeader(cluster)
    """
    Returns the updated size of cluster
    """

    re_size = r"[0-9]+$"

    num_sph = length(cluster["spheres"])

    return replace(
        cluster["line"],
        re_size => @sprintf "%2i" num_sph
        )
end

function UpdateCluster!(clusters, sph, rm_info)
    """
    Removes a given sphere from its cluster
    """

    cluster_id = sph.cluster_id

    sph_list = clusters[cluster_id]["spheres"]

    rm_idx = 0
    for (idx, s) in enumerate(sph_list)
        if s === sph
            rm_idx = idx
        end
    end

    deleteat!(clusters[cluster_id]["spheres"], rm_idx)
    rm_info[cluster_id] += 1
end

function FindNearest(c, all_spheres)

    d_arr = zeros(length(all_spheres))

    for (idx, s) in enumerate(all_spheres)
        d = sum((s.coordinates - c).^2)
        d_arr[idx] = d
    end

    return all_spheres[argmin(d_arr)]
end

function ParseSphere(line)
    vals = split(line)
    coords = [
        parse(Float64, c) for c in vals[2:4]
    ]
    return coords
end

function ParseMatchingSpheres(fn)

    re_cls = r"^cluster"
    re_sph = r"^ 9"

    headers = []
    clusters = Dict()

    current_cluster = 0
    sphere_id = 1

    all_spheres = []
    open(fn, "r") do io
        while !eof(io)
            line = readline(io)

            if occursin(re_cls, line)
                current_cluster += 1
                clusters[current_cluster] = Dict(
                    "line" => line,
                    "spheres" => []
                )

            elseif occursin(re_sph, line)
                coords = ParseSphere(line)
                sph = Sphere(sphere_id, current_cluster, coords, line)

                push!(all_spheres, sph)
                push!(clusters[current_cluster]["spheres"], sph)
                sphere_id += 1

            else
                push!(headers, line)
            end

        end
    end

    return headers, clusters, all_spheres
end

function ParseCoordinates(fn, rm_set)

    rm_coords = []
    open(fn, "r") do io
        while !eof(io)
            line = readline(io)
            ms_id = split(line)[end]
            if in(ms_id, rm_set)
                coords = ParseSphere(line)
                push!(rm_coords, coords)
            end
        end
    end
    return rm_coords
end

function main(args)

    ms_fn, coord_fn, rm_set = args

    headers, clusters, all_spheres = ParseMatchingSpheres(ms_fn)
    rm_coords = ParseCoordinates(coord_fn, rm_set)

    rm_info = zeros(Int32, length(clusters))
    for c in rm_coords
        closest_sph = FindNearest(c, all_spheres)
        UpdateCluster!(clusters, closest_sph, rm_info)
    end

    # print out headers
    for h in headers
        println(h)
    end

    # print out cluster spheres
    for idx in 1:length(clusters)

        cluster = clusters[idx]
        ch = GetClusterHeader(cluster)

        # cluster header
        println(ch)

        # sphere lines
        for sph in cluster["spheres"]
            println(sph.line)
        end

    end

    total_rm = sum(rm_info)
    println(stderr, "Removed : $total_rm")
    for (idx, rm_size) in enumerate(rm_info)
        println(stderr, "  From Cluster $idx : $rm_size")
    end

end

function print_help()
    println("Error : Expecting at least 3 positional arguments")
    println("julia PruneSpheres.jl <matching_spheres> <usage_coordinates> <indices> ... ")
end

function get_args()
    if length(ARGS) < 3
        print_help()
        exit()
    end

    matching_sphere_fn = ARGS[1]
    coordinate_usage = ARGS[2]
    rm_set = Set(ARGS[3:end])

    return (matching_sphere_fn, coordinate_usage, rm_set)
end

args = get_args()
main(args)

max_stars_repo_name: noamteyssier/ClusterDock

max_issues_repo_name: noamteyssier/ClusterDock

max_forks_repo_name: noamteyssier/ClusterDock


text:
# Autogenerated wrapper script for Libffi_jll for i686-w64-mingw32
export libffi

JLLWrappers.@generate_wrapper_header("Libffi")
JLLWrappers.@declare_library_product(libffi, "libffi-6.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libffi,
        "bin\\libffi-6.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Flux

# Zero Flux Given

using Base: tail
using Statistics, Random, LinearAlgebra
using Zygote, MacroTools, Juno, Reexport
using MacroTools: @forward
@reexport using NNlib
using Zygote: Params, @adjoint, gradient, pullback, @nograd

export gradient

export Chain, Dense, Maxout, RNN, LSTM, GRU, SamePad, Conv, CrossCor, ConvTranspose,
       AdaptiveMaxPool, AdaptiveMeanPool, GlobalMaxPool, GlobalMeanPool, MaxPool,
       MeanPool, flatten, DepthwiseConv, Dropout, AlphaDropout, LayerNorm, BatchNorm,
       InstanceNorm, GroupNorm, SkipConnection, params, fmap, cpu, gpu, f32, f64,
       testmode!, trainmode!

include("optimise/Optimise.jl")
using .Optimise
using .Optimise: @epochs
using .Optimise: skip
export Descent, ADAM, Momentum, Nesterov, RMSProp,
  ADAGrad, AdaMax, ADADelta, AMSGrad, NADAM, OADAM,
  ADAMW, RADAM, AdaBelief, InvDecay, ExpDecay,
  WeightDecay, ClipValue, ClipNorm

using CUDA
const use_cuda = Ref(false)

include("utils.jl")
include("zeros.jl")
include("onehot.jl")
include("functor.jl")

include("layers/stateless.jl")
include("layers/basic.jl")
include("layers/conv.jl")
include("layers/recurrent.jl")
include("layers/normalise.jl")

include("data/Data.jl")

include("losses/Losses.jl")
using .Losses # TODO: stop importing Losses in Flux's namespace in v0.12

include("deprecations.jl")

include("cuda/cuda.jl")

function __init__()
  use_cuda[] = CUDA.functional() # Can be overridden after load with `Flux.use_cuda[] = false`
  if CUDA.functional()
    if !CUDA.has_cudnn()
      @warn "CUDA.jl found cuda, but did not find libcudnn. Some functionality will not be available."
    end
  end
end

end # module


text:
@testset "ch2" begin
    Random.seed!(1234)
    rng = MersenneTwister(1234)
    function create_toy_data(func, sample_size, std)
        x = collect(range(0, stop = 1.0, length = sample_size))
        noise = rand(Normal(0.0, std), sample_size)
        return x, func(x) + noise
    end

    function sinusoidal(x)
        return sin.(2 * pi * x)
    end

    @testset "ch2_1" begin
        mu = [0.0, 1.0, 1.0, 1.0]
        bernoulli = BernoulliDist(mu)
        pdf(bernoulli, 0)
    end

    @testset "ch2_2" begin
        x = collect(range(0, stop = 1.0, length = 100))
        for (index, values) in enumerate([[0.1, 0.1], [1, 1], [2, 3], [8, 4]])
            a, b = values[1], values[2]
            beta = BetaDist(a, b)
            y = pdf(beta, x)
        end
    end

    @testset "ch2_3" begin
        x = collect(range(0, stop = 1.0, length = 100))

        beta = BetaDist(2, 2)

        bern = BernoulliDist(beta)
        fitting(bern, [1.0])
    end

    @testset "ch2_4" begin
        println("Maximum likehood estimation")

        model = BernoulliDist([1.0])
        count = draw(model, 10000)
        println("$(count) out of 10000 is label-1")

        println("Beysian estimation")
        model = BernoulliDist(BetaDist(1.0, 1.0))
        fitting(model, [1.0])
        count = draw(model, 10000)
        println("$(count) out of 10000 is label-1")
    end

    @testset "ch2_5" begin
        model = CategoricalDist([0.0])
        fitting(model, [[0.0 1.0 0.0]; [1.0 0.0 0.0]; [0.0 1.0 0.0]])
        println(model._mu)

        mu = DirichletDist(ones(3))
        model = CategoricalDist(mu)
        println("prior")
        println(model._dirichlet._alpha)
        trials = [[1.0 0.0 0.0]; [1.0 0.0 0.0]; [0.0 1.0 0.0]]
        fitting(model, trials)
        println("posterior")
        println(model._dirichlet._alpha)
    end

    @testset "ch2_6" begin
        dirichlet1 = DirichletDist([0.1, 0.1, 0.1])
        dirichlet2 = DirichletDist([1.0, 1.0, 1.0])
        dirichlet3 = DirichletDist([10.0, 10.0, 10.0])
    end

    @testset "ch2_7" begin
        mu_prior = GaussianDist(0.0, 0.1)
        model = GaussianBayesMeanDist(mu_prior, 0.1)

        x = collect((range(-1, stop = 1.0, length = 200)))

        fitting(model, [rand(Normal(0.8, 0.01))])

        fitting(model, [rand(Normal(0.8, 0.01))])
    end

    @testset "ch2_8" begin
        x = collect(range(0, stop = 2.0, length = 100))
        for (index, values) in enumerate([[0.1, 0.1], [1.0, 1.0], [2.0, 3.0], [4.0, 6.0]])
            gamma = GammaDist(values[1], values[2])
        end
    end

    @testset "ch2_9" begin
        X = rand(Uniform(1, 3), 2, 100)
        gaussian = MultivariateGaussianDist(2)
        fitting(gaussian, X)

        x = range(-5, stop = 10, length = 100)
        y = range(-5, stop = 10, length = 100)
        grids = [[i, j] for i in x, j in y]
        vals = [pdf(gaussian, grids[i, j]) for i = 1:100, j = 1:100]
    end
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
module NQCBase

using PeriodicTable
using Unitful, UnitfulAtomic
using StaticArrays
using Requires

include("unit_conversions.jl")
include("atoms.jl")
include("cells.jl")
include("io/extxyz.jl")

function __init__()
    @require PyCall="438e738f-606a-5dbb-bf0a-cddfbfd45ab0" @eval include("io/ase.jl")
end

end