text:
# these are just shorthads
using MonteCarlo: FVec64, FMat64, CVec64, CMat64

@testset "Field Cache & Interaction Matrix" begin
    T2 = NTuple{2}

    mc = DQMC(HubbardModel(8, 1, U = 1.0), field = DensityHirschField, beta = 1.0)
    @test mc.stack.field_cache isa MonteCarlo.StandardFieldCache{Float64, Float64, FVec64, Float64}
    @test MonteCarlo.interaction_eltype(mc.field) == Float64
    @test MonteCarlo.interaction_matrix_type(mc.field, mc.model) == Diagonal{Float64, FVec64}
    @test MonteCarlo.init_interaction_matrix(mc.field, mc.model) isa typeof(mc.stack.eV)

    mc = DQMC(HubbardModel(8, 1, U = -1.0), field = DensityHirschField, beta = 1.0)
    @test mc.stack.field_cache isa MonteCarlo.StandardFieldCache{ComplexF64, ComplexF64, CVec64, ComplexF64}
    @test MonteCarlo.interaction_eltype(mc.field) == ComplexF64
    @test MonteCarlo.interaction_matrix_type(mc.field, mc.model) == Diagonal{ComplexF64, CVec64}
    @test MonteCarlo.init_interaction_matrix(mc.field, mc.model) isa typeof(mc.stack.eV)

    mc = DQMC(HubbardModel(8, 1, U = 1.0), field = MagneticHirschField, beta = 1.0)
    @test mc.stack.field_cache isa MonteCarlo.StandardFieldCache{CVec64, CVec64, T2{CVec64}, ComplexF64}
    @test MonteCarlo.interaction_eltype(mc.field) == ComplexF64
    @test MonteCarlo.interaction_matrix_type(mc.field, mc.model) == Diagonal{ComplexF64, CVec64}
    @test MonteCarlo.init_interaction_matrix(mc.field, mc.model) isa typeof(mc.stack.eV)

    mc = DQMC(HubbardModel(8, 1, U = -1.0), field = MagneticHirschField, beta = 1.0)
    @test mc.stack.field_cache isa MonteCarlo.StandardFieldCache{FVec64, FVec64, T2{FVec64}, Float64}
    @test MonteCarlo.interaction_eltype(mc.field) == Float64
    @test MonteCarlo.interaction_matrix_type(mc.field, mc.model) == Diagonal{Float64, FVec64}
    @test MonteCarlo.init_interaction_matrix(mc.field, mc.model) isa typeof(mc.stack.eV)

    mc = DQMC(HubbardModel(8, 1, U = 1.0), field = MagneticGHQField, beta = 1.0)
    @test mc.stack.field_cache isa MonteCarlo.StandardFieldCache{CVec64, CVec64, T2{CVec64}, ComplexF64}
    @test MonteCarlo.interaction_eltype(mc.field) == ComplexF64
    @test MonteCarlo.interaction_matrix_type(mc.field, mc.model) == Diagonal{ComplexF64, CVec64}
    @test MonteCarlo.init_interaction_matrix(mc.field, mc.model) isa typeof(mc.stack.eV)

    mc = DQMC(HubbardModel(8, 1, U = -1.0), field = MagneticGHQField, beta = 1.0)
    @test mc.stack.field_cache isa MonteCarlo.StandardFieldCache{FVec64, FVec64, T2{FVec64}, Float64}
    @test MonteCarlo.interaction_eltype(mc.field) == Float64
    @test MonteCarlo.interaction_matrix_type(mc.field, mc.model) == Diagonal{Float64, FVec64}
    @test MonteCarlo.init_interaction_matrix(mc.field, mc.model) isa typeof(mc.stack.eV)
end

@testset "Lookup tables" begin
    param = MonteCarlo.DQMCParameters(beta = 1.0)
    model = HubbardModel()
    field = MagneticGHQField(param, model)

    @test field.α == sqrt(-0.5 * 0.1 * ComplexF64(model.U))

    # See ALF Documentation for the formulas of η and γ
    @test field.γ[1] ≈ Float64(BigFloat(1) - sqrt(BigFloat(6)) / BigFloat(3))  rtol=1e-15
    @test field.γ[2] ≈ Float64(BigFloat(1) + sqrt(BigFloat(6)) / BigFloat(3))  rtol=1e-15
    @test field.γ[3] ≈ Float64(BigFloat(1) + sqrt(BigFloat(6)) / BigFloat(3))  rtol=1e-15
    @test field.γ[4] ≈ Float64(BigFloat(1) - sqrt(BigFloat(6)) / BigFloat(3))  rtol=1e-15

    @test field.η[1] ≈ Float64(- sqrt(BigFloat(2) * (BigFloat(3) + sqrt(BigFloat(6)))))  rtol=1e-15
    @test field.η[2] ≈ Float64(- sqrt(BigFloat(2) * (BigFloat(3) - sqrt(BigFloat(6)))))  rtol=1e-15
    @test field.η[3] ≈ Float64(+ sqrt(BigFloat(2) * (BigFloat(3) - sqrt(BigFloat(6)))))  rtol=1e-15
    @test field.η[4] ≈ Float64(+ sqrt(BigFloat(2) * (BigFloat(3) + sqrt(BigFloat(6)))))  rtol=1e-15

    @test field.choices[1, :] == Int8[2, 3, 4]
    @test field.choices[2, :] == Int8[1, 3, 4]
    @test field.choices[3, :] == Int8[1, 2, 4]
    @test field.choices[4, :] == Int8[1, 2, 3]
end

@testset "Interface" begin
    fields = (MagneticHirschField, DensityHirschField, MagneticGHQField)

    for field_type in fields
        mc = DQMC(HubbardModel(2, 2), beta = 1.0, field = field_type)
        f = mc.field

        @test MonteCarlo.maybe_to_float(ComplexF64(1, 0)) === 1.0
        @test MonteCarlo.maybe_to_float(ComplexF64(0, 1)) === ComplexF64(0, 1)

        @test MonteCarlo.field(mc) == f
        @test MonteCarlo.conf(f) == f.conf
        @test MonteCarlo.temp_conf(f) == f.temp_conf
        @test MonteCarlo.length(f) == length(f.conf)

        c = rand(f)
        @test MonteCarlo.conf(f) != c
        MonteCarlo.conf!(f, c)
        @test MonteCarlo.conf(f) == c
        rand!(f)
        @test MonteCarlo.conf(f) != c

        compressed = MonteCarlo.compress(f)
        @test compressed isa BitArray
        @test MonteCarlo.compressed_conf_type(f) <: BitArray
        @test MonteCarlo.decompress(f, compressed) == MonteCarlo.conf(f)
        c = MonteCarlo.conf(f); rand!(f)
        MonteCarlo.decompress!(f, compressed)
        @test MonteCarlo.conf(f) == c

        @test MonteCarlo.nflavors(f) == (2 - (f isa DensityHirschField))
        E = f isa DensityHirschField ? f.α * sum(f.conf) : 0.0
        @test MonteCarlo.energy_boson(f) == E
    end
end

using StructArrays, LinearAlgebra
using MonteCarlo, Test

@testset "Linear Algebra" begin
    i = 3; N = 4; C = ComplexF64; SA = StructArray

    G = rand(N, N)
    cache = MonteCarlo.StandardFieldCache(rand(), rand(), rand(), rand(N), rand(N), rand())
    cache.R = 1 + cache.Δ * (1 - G[i, i])
    cache.detratio = cache.R * cache.R

    invRΔ = cache.Δ / cache.R
    MonteCarlo.vldiv22!(cache, cache.R, cache.Δ)
    @test invRΔ ≈ cache.invRΔ rtol = 1e-14

    IG = (I - Matrix(G))[:, i:N:end]
    MonteCarlo.vsub!(cache.IG, I, G, i, N)
    # @test cache.IG ≈ IG rtol = 1e-14

    g = invRΔ * Matrix(G)[i:N:end, :]
    MonteCarlo.vmul!(cache.G, cache.invRΔ, G, i, N)
    # @test g ≈ cache.G rtol = 1e-14

    Q = Matrix(G) - IG * g
    MonteCarlo.vsubkron!(G, cache.IG, cache.G)
    @test G ≈ Q rtol = 1e-14

    G = SA(rand(ComplexF64, N, N))
    cache = MonteCarlo.StandardFieldCache(rand(C), rand(C), rand(C), SA(rand(C, N)), SA(rand(C, N)), rand(C))
    cache.R = 1 + cache.Δ * (1 - G[i, i])
    cache.detratio = cache.R * cache.R

    invRΔ = cache.Δ / cache.R
    MonteCarlo.vldiv22!(cache, cache.R, cache.Δ)
    @test invRΔ ≈ cache.invRΔ rtol = 1e-14

    IG = (I - Matrix(G))[:, i:N:end]
    MonteCarlo.vsub!(cache.IG, I, G, i, N)
    # @test cache.IG ≈ IG rtol = 1e-14

    g = invRΔ * Matrix(G)[i:N:end, :]
    MonteCarlo.vmul!(cache.G, cache.invRΔ, G, i, N)
    # @test g ≈ cache.G rtol = 1e-14

    Q = Matrix(G) - IG * g
    MonteCarlo.vsubkron!(G, cache.IG, cache.G)
    @test G ≈ Q rtol = 1e-14

    G = rand(2N, 2N)
    cache = MonteCarlo.StandardFieldCache(rand(2), rand(2, 2), rand(2, 2), rand(2N, 2), rand(2N, 2), rand())
    cache.R .= I + Diagonal(cache.Δ) * (I - G[i:N:end, i:N:end])
    cache.detratio = det(cache.R)

    invRΔ = inv(cache.R) * Diagonal(cache.Δ)
    MonteCarlo.vldiv22!(cache, cache.R, cache.Δ)
    @test invRΔ ≈ cache.invRΔ rtol = 1e-14

    IG = (I - Matrix(G))[:, i:N:end]
    MonteCarlo.vsub!(cache.IG, I, G, i, N)
    # @test cache.IG ≈ IG rtol = 1e-14

    g = invRΔ * Matrix(G)[i:N:end, :]
    MonteCarlo.vmul!(cache.G, cache.invRΔ, G, i, N)
    # @test g ≈ cache.G rtol = 1e-14

    Q = Matrix(G) - IG * g
    MonteCarlo.vsubkron!(G, cache.IG, cache.G)
    @test G ≈ Q rtol = 1e-14

    G = MonteCarlo.BlockDiagonal(rand(N, N), rand(N, N))
    cache = MonteCarlo.StandardFieldCache(rand(2), rand(2), rand(2), (rand(N), rand(N)), (rand(N), rand(N)), rand())
    cache.R[1] = 1 + cache.Δ[1] * (1 - G[i, i])
    cache.R[2] = 1 + cache.Δ[2] * (1 - G[i+N, i+N])
    cache.detratio = cache.R[1] * cache.R[2]

    invRΔ = cache.Δ ./ cache.R
    MonteCarlo.vldiv22!(cache, cache.R, cache.Δ)
    @test invRΔ ≈ cache.invRΔ rtol = 1e-14

    IG = (I - Matrix(G))[:, i:N:end]
    MonteCarlo.vsub!(cache.IG, I, G, i, N)
    # @test cache.IG ≈ IG rtol = 1e-14

    g = Diagonal(invRΔ) * Matrix(G)[i:N:end, :]
    MonteCarlo.vmul!(cache.G, cache.invRΔ, G, i, N)
    # @test g ≈ cache.G rtol = 1e-14

    Q = Matrix(G) - IG * g
    MonteCarlo.vsubkron!(G, cache.IG, cache.G)
    @test G ≈ Q rtol = 1e-14

    G = MonteCarlo.BlockDiagonal(SA(rand(ComplexF64, N, N)), SA(rand(ComplexF64, N, N)))
    cache = MonteCarlo.StandardFieldCache(SA(rand(C, 2)), SA(rand(C, 2)), SA(rand(C, 2)), (SA(rand(C, N)), SA(rand(C, N))), (SA(rand(C, N)), SA(rand(C, N))), rand(C))
    cache.R[1] = 1 + cache.Δ[1] * (1 - G[i, i])
    cache.R[2] = 1 + cache.Δ[2] * (1 - G[i+N, i+N])
    cache.detratio = cache.R[1] * cache.R[2]

    invRΔ = cache.Δ ./ cache.R
    MonteCarlo.vldiv22!(cache, cache.R, cache.Δ)
    @test invRΔ ≈ cache.invRΔ rtol = 1e-14

    IG = (I - Matrix(G))[:, i:N:end]
    MonteCarlo.vsub!(cache.IG, I, G, i, N)
    # @test cache.IG ≈ IG rtol = 1e-14

    g = Diagonal(invRΔ) * Matrix(G)[i:N:end, :]
    MonteCarlo.vmul!(cache.G, cache.invRΔ, G, i, N)
    # @test g ≈ cache.G rtol = 1e-14

    Q = Matrix(G) - IG * g
    MonteCarlo.vsubkron!(G, cache.IG, cache.G)
    @test G ≈ Q rtol = 1e-14
end

@testset "Sign Problem in field - model Combinations" begin
    models = (HubbardModel(8, 1, U = 1.0), HubbardModel(8, 1, U = -1.0))
    fields = (MagneticHirschField, DensityHirschField, MagneticGHQField)

    for m in models, field in fields
        mc = DQMC(m, field = field, beta=3.1, thermalization=50, sweeps = 50, print_rate = 100)
        run!(mc, verbose=false)

        @testset "$(nameof(field)) $(nameof(typeof(m)))" begin
            @test mc.analysis.imaginary_probability.count == 0
            @test mc.analysis.negative_probability.count == 0
            @test mc.analysis.propagation_error.count == 0
        end
    end
end


text:
#=
  Decomposition of global constraint global_cardinality in Julia ConstraintSolver.

  See Global Constraint Catalog
  http://www.emn.fr/x-info/sdemasse/gccat/Cglobal_cardinality.html
  """
  Example:
  <3,3,8,6>,
   (val-3 noccurrence-2,
    val−5 noccurrence-0,
    val-6 noccurrence-1)

  The global_cardinality constraint holds since values 3, 5 and 6
  respectively occur 2, 0 and 1 times within the collection
  <3,3,8,6> and since no constraint was specified for value 8.
  """

  Notes:
     - The version here is limited. See below for details.
     - Picat has already a global_cardinality/2, but the details
       of the parameters are different.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK
const CS = ConstraintSolver
include("constraints_utils.jl")

#
# Note: global_cardinality_count is defined in constraints_utils.jl
#
function gcc_test(n=5,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS, # <-
                                                            # "traverse_strategy"=>:DBFS,

                                                            "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            # "branch_split"=>:InHalf, # <-

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>16,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                        ))

    @variable(model, 1 <= a[1:n] <= n, Int)
    @variable(model, 0 <= gcc[1:n] <= n, Int)

    global_cardinality_count(model, a, gcc)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                ax = convert.(Integer,JuMP.value.(a; result=sol))
                gccx = convert.(Integer,JuMP.value.(gcc; result=sol))
                # println("position:$positionx\n")
                println("a:$ax gcc:$gccx\n")
            end
        end
    else
        println("status:$status")
    end
    return status
end

@time gcc_test(5)

ext: jl

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
export MonoidalCategoryAdditive, SymmetricMonoidalCategoryAdditive,
  FreeSymmetricMonoidalCategoryAdditive, oplus, ⊕, mzero, swap,
  MonoidalCategoryWithCodiagonals, CocartesianCategory, FreeCocartesianCategory,
  plus, zero, copair, coproj1, coproj2,
  MonoidalCategoryWithBidiagonalsAdditive, SemiadditiveCategory, AdditiveCategory,
  mcopy, delete, pair, proj1, proj2, Δ, ◊, +, antipode,
  HypergraphCategoryAdditive

import Base: collect, ndims, +, zero

# Monoidal category
###################

""" Theory of *monoidal categories*, in additive notation

Mathematically the same as [`MonoidalCategory`](@ref) but with different
notation.
"""
@signature MonoidalCategoryAdditive{Ob,Hom} <: Category{Ob,Hom} begin
  oplus(A::Ob, B::Ob)::Ob
  oplus(f::(A → B), g::(C → D))::((A ⊕ C) → (B ⊕ D)) <=
    (A::Ob, B::Ob, C::Ob, D::Ob)
  @op (⊕) := oplus
  mzero()::Ob
end

# Convenience constructors
oplus(xs::Vector{T}) where T = isempty(xs) ? mzero(T) : foldl(oplus, xs)
oplus(x, y, z, xs...) = oplus([x, y, z, xs...])

# Overload `collect` and `ndims` as for multiplicative monoidal categories.
collect(expr::ObExpr{:oplus}) = vcat(map(collect, args(expr))...)
collect(expr::ObExpr{:mzero}) = roottypeof(expr)[]
ndims(expr::ObExpr{:oplus}) = sum(map(ndims, args(expr)))
ndims(expr::ObExpr{:mzero}) = 0

function show_unicode(io::IO, expr::Union{ObExpr{:oplus},HomExpr{:oplus}}; kw...)
  Syntax.show_unicode_infix(io, expr, "⊕"; kw...)
end
show_unicode(io::IO, expr::ObExpr{:mzero}; kw...) = print(io, "O")

function show_latex(io::IO, expr::Union{ObExpr{:oplus},HomExpr{:oplus}}; kw...)
  Syntax.show_latex_infix(io, expr, "\\oplus"; kw...)
end
show_latex(io::IO, expr::ObExpr{:mzero}; kw...) = print(io, "O")

# Symmetric monoidal category
#############################

""" Theory of *symmetric monoidal categories*, in additive notation

Mathematically the same as [`SymmetricMonoidalCategory`](@ref) but with
different notation.
"""
@signature SymmetricMonoidalCategoryAdditive{Ob,Hom} <:
    MonoidalCategoryAdditive{Ob,Hom} begin
  swap(A::Ob, B::Ob)::Hom(oplus(A,B),oplus(B,A))
end

@syntax FreeSymmetricMonoidalCategoryAdditive{ObExpr,HomExpr} SymmetricMonoidalCategoryAdditive begin
  oplus(A::Ob, B::Ob) = associate_unit(new(A,B), mzero)
  oplus(f::Hom, g::Hom) = associate(new(f,g))
  compose(f::Hom, g::Hom) = associate(new(f,g; strict=true))
end

function show_latex(io::IO, expr::HomExpr{:swap}; kw...)
  Syntax.show_latex_script(io, expr, "\\sigma")
end

# Cocartesian category
######################

""" Theory of *monoidal categories with codiagonals*

A monoidal category with codiagonals is a symmetric monoidal category equipped
with coherent collections of merging and creating morphisms (monoids).
Unlike in a cocartesian category, the naturality axioms need not be satisfied.

For references, see [`MonoidalCategoryWithDiagonals`](@ref).
"""
@theory MonoidalCategoryWithCodiagonals{Ob,Hom} <:
    SymmetricMonoidalCategoryAdditive{Ob,Hom} begin
  plus(A::Ob)::((A ⊕ A) → A)
  zero(A::Ob)::(mzero() → A)

  # Commutative monoid axioms.
  (plus(A) ⊕ id(A)) ⋅ plus(A) == (id(A) ⊕ plus(A)) ⋅ plus(A) ⊣ (A::Ob)
  (zero(A) ⊕ id(A)) ⋅ plus(A) == id(A) ⊣ (A::Ob)
  (id(A) ⊕ zero(A)) ⋅ plus(A) == id(A) ⊣ (A::Ob)
  plus(A) == swap(A,A) ⋅ plus(A) ⊣ (A::Ob)

  # Coherence axioms.
  plus(A⊕B) == (id(A) ⊕ swap(B,A) ⊕ id(B)) ⋅ (plus(A) ⊕ plus(B)) ⊣ (A::Ob, B::Ob)
  zero(A⊕B) == zero(A) ⊕ zero(B) ⊣ (A::Ob, B::Ob)
  plus(mzero()) == id(mzero())
  zero(mzero()) == id(mzero())
end

""" Theory of *cocartesian (monoidal) categories*

For the traditional axiomatization of coproducts, see
[`CategoryWithCoproducts`](@ref).
"""
@theory CocartesianCategory{Ob,Hom} <: MonoidalCategoryWithCodiagonals{Ob,Hom} begin
  copair(f::(A → C), g::(B → C))::((A ⊕ B) → C) <= (A::Ob, B::Ob, C::Ob)
  coproj1(A::Ob, B::Ob)::(A → (A ⊕ B))
  coproj2(A::Ob, B::Ob)::(B → (A ⊕ B))

  # Definitions of copairing and coprojections.
  copair(f,g) == (f⊕g)⋅plus(C) ⊣ (A::Ob, B::Ob, C::Ob, f::(A → C), g::(B → C))
  coproj1(A,B) == id(A)⊕zero(B) ⊣ (A::Ob, B::Ob)
  coproj2(A,B) == zero(A)⊕id(B) ⊣ (A::Ob, B::Ob)

  # Naturality axioms.
  plus(A)⋅f == (f⊕f)⋅plus(B) ⊣ (A::Ob, B::Ob, f::(A → B))
  zero(A)⋅f == zero(B) ⊣ (A::Ob, B::Ob, f::(A → B))
end

""" Syntax for a free cocartesian category.

In this syntax, the copairing and inclusion operations are defined using merging
and creation, and do not have their own syntactic elements. This convention
could be dropped or reversed.
"""
@syntax FreeCocartesianCategory{ObExpr,HomExpr} CocartesianCategory begin
  oplus(A::Ob, B::Ob) = associate_unit(new(A,B), mzero)
  oplus(f::Hom, g::Hom) = associate(new(f,g))
  compose(f::Hom, g::Hom) = associate(new(f,g; strict=true))

  copair(f::Hom, g::Hom) = compose(oplus(f,g), plus(codom(f)))
  coproj1(A::Ob, B::Ob) = oplus(id(A), zero(B))
  coproj2(A::Ob, B::Ob) = oplus(zero(A), id(B))
end

function show_latex(io::IO, expr::HomExpr{:plus}; kw...)
  if length(args(expr)) >= 2
    Syntax.show_latex_infix(io, expr, "+"; kw...)
  else
    Syntax.show_latex_script(io, expr, "\\nabla")
  end
end

function show_latex(io::IO, expr::HomExpr{:zero}; kw...)
  Syntax.show_latex_script(io, expr, "0")
end

# Additive category
###################

""" Theory of *monoidal categories with bidiagonals*, in additive notation

Mathematically the same as [`MonoidalCategoryWithBidiagonals`](@ref) but written
additively, instead of multiplicatively.
"""
@theory MonoidalCategoryWithBidiagonalsAdditive{Ob,Hom} <:
    MonoidalCategoryWithCodiagonals{Ob,Hom} begin
  mcopy(A::Ob)::(A → (A ⊕ A))
  @op (Δ) := mcopy
  delete(A::Ob)::(A → mzero())
  @op (◊) := delete

  # Commutative comonoid axioms.
  Δ(A) == Δ(A) ⋅ swap(A,A) ⊣ (A::Ob)
  Δ(A) ⋅ (Δ(A) ⊕ id(A)) == Δ(A) ⋅ (id(A) ⊕ Δ(A)) ⊣ (A::Ob)
  Δ(A) ⋅ (◊(A) ⊕ id(A)) == id(A) ⊣ (A::Ob)
  Δ(A) ⋅ (id(A) ⊕ ◊(A)) == id(A) ⊣ (A::Ob)
end

""" Theory of *semiadditive categories*

Mathematically the same as [`BiproductCategory`](@ref) but written additively,
instead of multiplicatively.
"""
@theory SemiadditiveCategory{Ob,Hom} <:
    MonoidalCategoryWithBidiagonalsAdditive{Ob,Hom} begin
  pair(f::(A → B), g::(A → C))::(A → (B ⊕ C)) ⊣ (A::Ob, B::Ob, C::Ob)
  copair(f::(A → C), g::(B → C))::((A ⊕ B) → C) ⊣ (A::Ob, B::Ob, C::Ob)
  proj1(A::Ob, B::Ob)::((A ⊕ B) → A)
  proj2(A::Ob, B::Ob)::((A ⊕ B) → B)
  coproj1(A::Ob, B::Ob)::(A → (A ⊕ B))
  coproj2(A::Ob, B::Ob)::(B → (A ⊕ B))

  plus(f::(A → B), g::(A → B))::(A → B) ⊣ (A::Ob, B::Ob)
  @op (+) := plus

  # Naturality axioms.
  f⋅Δ(B) == Δ(A)⋅(f⊕f) ⊣ (A::Ob, B::Ob, f::(A → B))
  f⋅◊(B) == ◊(A) ⊣ (A::Ob, B::Ob, f::(A → B))
  plus(A)⋅f == (f⊕f)⋅plus(B) ⊣ (A::Ob, B::Ob, f::(A → B))
  zero(A)⋅f == zero(B) ⊣ (A::Ob, B::Ob, f::(A → B))

  # Bimonoid axioms. (These follow from naturality + coherence axioms.)
  plus(A)⋅Δ(A) == (Δ(A)⊕Δ(A)) ⋅ (id(A)⊕swap(A,A)⊕id(A)) ⋅ (plus(A)⊕plus(A)) ⊣ (A::Ob)
  plus(A)⋅◊(A) == ◊(A) ⊕ ◊(A) ⊣ (A::Ob)
  zero(A)⋅Δ(A) == zero(A) ⊕ zero(A) ⊣ (A::Ob)
  zero(A)⋅◊(A) == id(mzero()) ⊣ (A::Ob)
end

""" Theory of *additive categories*

An additive category is a biproduct category enriched in abelian groups. Thus,
it is a semiadditive category where the hom-monoids have negatives.
"""
@theory AdditiveCategory{Ob,Hom} <: SemiadditiveCategory{Ob,Hom} begin
  antipode(A::Ob)::(A → A)

  # Antipode axioms.
  antipode(A) ⋅ f == f ⋅ antipode(B) ⊣ (A::Ob, B::Ob, f::(A → B))
  Δ(A)⋅(id(A)⊕antipode(A))⋅plus(A) == ◊(A)⋅zero(A) ⊣ (A::Ob)
  Δ(A)⋅(antipode(A)⊕id(A))⋅plus(A) == ◊(A)⋅zero(A) ⊣ (A::Ob)
end

# Hypergraph category
#####################

""" Theory of *hypergraph categories*, in additive notation

Mathematically the same as [`HypergraphCategory`](@ref) but with different
notation.
"""
@signature HypergraphCategoryAdditive{Ob,Hom} <:
    SymmetricMonoidalCategoryAdditive{Ob,Hom} begin
  # Supply of Frobenius monoids.
  mcopy(A::Ob)::(A → (A ⊕ A))
  @op (Δ) := mcopy
  delete(A::Ob)::(A → mzero())
  @op (◊) := delete
  mmerge(A::Ob)::((A ⊕ A) → A)
  @op (∇) := mmerge
  create(A::Ob)::(mzero() → A)
  @op (□) := create

  # Self-dual compact closed category.
  dunit(A::Ob)::(mzero() → (A ⊕ A))
  dcounit(A::Ob)::((A ⊕ A) → mzero())
  dagger(f::(A → B))::(B → A) ⊣ (A::Ob, B::Ob)
end


text:
#Written By Johnathan Bizzano
export AscendingSorter, DescendingSorter, DefSorter, RandomSorter

"""
    Sorting Library for the Matrix DataSet Object
"""

"""
    Sort the data set from least to greatest
"""
function AscendingSorter(matrix::MatrixDataSet)

end

"""
    Sort the data set from greatest to least
"""
function DescendingSorter(matrix::MatrixDataSet)

end

"""
    Sort Randomly
"""
function RandomSorter(matrix::MatrixDataSet)

end


text:
module LinearAlgebra

using Reexport

include("GLA.jl")
include("StructuredGLA.jl")
include("bayes.jl")

@reexport using .GraphicalLinearAlgebra
@reexport using .StructuredGraphicalLinearAlgebra
@reexport using .Bayes
end


text:
using VITA2D

L = 12
P = 1
θ = [0.2, 0.2, 0.0]
# θ = rand(3P)

sgd_params = VITA2D.SGDParams(lr=0.2, N_mc=10^3, N_sgd=200, num_avg=10)

E = sgd(L=L,
        P=P,
        hX = 3.0,
        θ = θ,
        SGD=sgd_params,
        print=true)


text:
formatcode(::Type{format"WAV"}) = SF_FORMAT_WAV
formatcode(::Type{format"FLAC"}) = SF_FORMAT_FLAC
formatcode(::Type{format"OGG"}) = SF_FORMAT_OGG

subformatcode(::Type{PCM16Sample}) = SF_FORMAT_PCM_16
subformatcode(::Type{PCM32Sample}) = SF_FORMAT_PCM_32
subformatcode(::Type{Float32}) = SF_FORMAT_FLOAT
subformatcode(::Type{Float64}) = SF_FORMAT_DOUBLE

"""Take a LibSndFile format code and return a suitable sample type"""
function fmt_to_type(fmt)
  mapping = Dict{UInt32, Type}(
                               SF_FORMAT_PCM_S8 => PCM16Sample,
                               SF_FORMAT_PCM_U8 => PCM16Sample,
                               SF_FORMAT_PCM_16 => PCM16Sample,
                               SF_FORMAT_PCM_24 => PCM32Sample,
                               SF_FORMAT_PCM_32 => PCM32Sample,
                               SF_FORMAT_FLOAT => Float32,
                               SF_FORMAT_DOUBLE => Float64,
                               SF_FORMAT_VORBIS => Float32,
                              )

  masked = fmt & SF_FORMAT_SUBMASK
  masked in keys(mapping) || error("Format code $masked not recognized by LibSndFile.jl")

  mapping[masked]
end

mutable struct SF_INFO
  frames::sf_count_t
  samplerate::Int32
  channels::Int32
  format::Int32
  sections::Int32
  seekable::Int32
end

SF_INFO() = SF_INFO(0, 0, 0, 0, 0, 0)

function sf_open(fname::String, mode, sfinfo)
  ## this fixes #34 however is unstable and breaks test randomly
  ## it's difficult to debug it without a Windows machine
  #ptr = pointer(transcode(Cwchar_t,fname))
  #filePtr = ccall((:sf_wchar_open, libsndfile), Ptr{Cvoid},
  #                (Cwstring, Int32, Ref{SF_INFO}),
  #                Cwstring(ptr), mode, sfinfo
  filePtr = ccall((:sf_open, libsndfile), Ptr{Cvoid},
                  (Cstring, Int32, Ref{SF_INFO}),
                  fname, mode, sfinfo)
  if filePtr == C_NULL
    error("LibSndFile.jl error while opening $fname: ", sf_strerror(C_NULL))
  end

  filePtr
end

function sf_close(filePtr)
  err = ccall((:sf_close, libsndfile), Int32, (Ptr{Cvoid},), filePtr)
  if err != 0
    error("LibSndFile.jl error: Failed to close file: ", sf_strerror(filePtr))
  end
end

sf_seek(filePtr, frames::sf_count_t, whence::Integer) =
ccall((:sf_seek, libsndfile), Int64,
      (Ptr{Cvoid}, Int64, Int32),
      filePtr, frames, whence)


text:
using SQLdf
using Test
import DataFrames: DataFrame, rename

@testset "SQLdf.jl" begin
    Main.eval(Meta.parse(""" x = DataFrame((a=1:14,  b=14:-1:1, c =split("Julia is great",""))) """))

    qr = sqldf("""
               select *
               from x
               where a <= 3
               order by a
              """)
    @test size(qr) == (3,3)
    @test qr[1,1] == 1

    qr = sqldf("""
               select *
               from x
               where a <= 3 and c = "J"
               order by a
              """)

    @test size(qr) == (1,3)
    @test qr[1,2] == 14

    Main.eval(Meta.parse(""" y = DataFrame((a=1:14, c=split("Julia is fast!",""))) """))

    qr = sqldf("""
               select x.*, y.*
               from x join y on x.c = y.c
               order by x.a
              """)

    qr = sqldf("""
               select *
               from x join y on x.c = y.c
               order by a
              """)
    @test names(qr) == ["a", "b", "c", "a:1", "c:1"]
    @test qr[1,3] == "J"
    @test qr[1,5] == "J"
    @test qr[18,5] == "t"

    q = "select * from x"
    @test sqldf(q) == @sqldf q
    @test sqldf("select * from x") == @sqldf "select * from x"

end


text:
module WCvsLHE

using Images, ImageFiltering
using FFTW
import OffsetArrays

const Lift{T} = Array{T,3}
const Kern{T,n} = OffsetArrays.OffsetArray{T,n,Array{T,n}}

include("lift.jl")

export wc, lhe, Params

normalize(x::Array) = typeof(x)(x - ones(size(x))*minimum(x))/(maximum(x)-minimum(x))

function project(F::Lift; normalized::Bool = true, args...)
    x = dropdims(sum(F, dims = 3), dims = 3)
    normalized ? x = normalize(x) : nothing
    return x
end

# Approximation of sigmoid for LHE
# Coefficients for the polynomial approximation of the sigmoid
const c = reverse([-7.7456e+00,
                3.1255e-16,
                1.5836e+01,
                -1.8371e-15,
                -1.1013e+01,
                4.4531e-16,
                3.7891e+00,
                1.2391e-15]);

#Degree of the approximating polynomial
const n = length(c)-1;

gr(i::Int, img::Array, ω::Kern) = imfilter(img.^i,ω)
a(i::Int, img::Array)=sum( [(-1.)^(j-i+1)*c[j+1]*binomial(j,i)*img.^(j-i) for j in i:n] )
R(img::Array, ω::Kern) = sum( [ a(i,img) .* gr(i,img,ω) for i in 0:n ] )

# Weight functions generation
function adapted_gaussian(σ::Real, I::Array)
    # Generates a Gaussian with variance σ and correct
    # dimensions w.r.t. to I
    if ndims(I) == 2
        var = (σ, σ)
    else
        var = (σ, σ, (size(I, 3)/size(I, 1))*σ)
    end

    # var = tuple([σ for _ in 1:ndims(I)]...)
    Kernel.gaussian(var) |> reflect # Convolution not correlation
end

# Local Mean Average
LMA(σμ::Real, I0::Array{T, 2}) where T = imfilter(I0, adapted_gaussian(σμ, I0))

# Type to save the results
struct Result{n}
    res :: Array{Float64, n}
    iter :: Int
    tol :: Float64
end
Base.getindex(R::Result, etc...)  = getindex(R.res, etc...)
iterations(R::Result) = R.iter
tolerance(R::Result) = R.tol
project(R::Result{3}, etc...) = Result{2}(project(R.res, etc...), iterations(R), tolerance(R))
show(R::Result{2}) = Gray.(R[:,:])

# Actual procedure
function gradient_descent(I0::Array, W, λ::Real, lma::Array; Δt::Real = .15, threshold::Real = .01, max_iter::Int = 50, M::Real = 1, args...)
    # threshold : the iterations stops when two successive terms have a weighted L2 difference smaller than this quantity
    # max_iter: maximum number of iterations
    next(Ik) = (1-(1+λ)Δt)Ik + (lma + W(Ik)/(2*M))Δt + λ*Δt*I0
    proj(x) = sum(x, dims = 3) # it does nothing if x has 2 dims
    diff(a,b) = sqrt( sum(( proj(a-b) ).^2)/sum(proj(a).^2) )

    prec = I0
    cur = next(I0)
    i = 1
    while diff(prec, cur) >= threshold && i < max_iter
        (prec, cur) = (cur, next(cur))
        i += 1
    end

    Result{ndims(I0)}(cur, i, diff(prec, cur))
end

function algo(I0::Array, σμ::Real, σw::Real, λ::Real; algo_type = :planar, model = :wc, θs::Int = 30, α::Int = 5, verbose::Bool = false, args...)

	if algo_type == :cortical
		F0 = lift(I0, θs; args...)
		lma = lift(LMA(σμ, I0), θs; args...)
	elseif algo_type == :planar
		F0 = copy(I0)
		lma = LMA(σμ, I0)
	else
		error("Type of evolution must be :cortical or :planar")
	end

	if model == :wc
        function σ(x)
            if x < 1/2 - 1/α
                y = 1
            elseif x > 1/2 + 1/α
                y = -1
            else
                y = (1/2-x)*α
            end
            y
        end

		W(I) = imfilter(σ.(I), adapted_gaussian(σw, F0))
        res = gradient_descent(F0, W, λ, lma; args...)
	elseif model == :lhe
		W2(I) = R(I, adapted_gaussian(σw, F0))
        res = gradient_descent(F0, W2, λ, lma; args...)
    else
		error("Model must be :wc or :lhe")
	end

	if algo_type == :cortical
		res = project(res)
	end

    verbose ? res : show(res)
end

wc(I0, σμ, σw, λ; args...) = algo(I0, σμ, σw, λ, model = :wc; args...)
lhe(I0, σμ, σw, λ; args...) = algo(I0, σμ, σw, λ, model = :lhe; args...)

# A useful shorthand:
struct Params
    σμ :: Float64 # std deviation of local mean average
    σw :: Float64 # std deviation of interaction kernel
    λ :: Float64 # attachment to the data
    M :: Float64 # the interaction weight is ν=1/(2M)
end
σμ(p::Params) = p.σμ
σw(p::Params) = p.σw
λ(p::Params) = p.λ
M(p::Params) = p.M

wc(I, p::Params; args...) = wc(Float64.(I), σμ(p), σw(p), λ(p), M = M(p); args...)
lhe(I, p::Params; args...) = lhe(Float64.(I), σμ(p), σw(p), λ(p), M = M(p); args...)

end

max_stars_repo_name: dprn/WCvsLHE

max_issues_repo_name: dprn/WCvsLHE

max_forks_repo_name: dprn/WCvsLHE


text:
using IScatterSpectrum
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# Scattering volume for 230 MHz monostatic along the magnetic field
sv = IScatterSpectrum.ScatterVolume(230e6, 1e-6, 50000e-9, 0.0)

# Plasma, typical Ne, Te, and Ti
p  = IScatterSpectrum.Plasma(1.5e11, 2000e0, 1000e0, sv)

# Frequency array good for 230 MHz
f = 5.0:5:5000

s = [IScatterSpectrum.pwrspec(p, sv); IScatterSpectrum.pwrspec(1.0, p, sv)]
@test abs(s[2]-s[1])<1e-5


text:
function always_assert(cond::Bool, msg::AbstractString)::Nothing
    _msg::String = convert(String, msg)::String
    cond || throw(AlwaysAssertionError(_msg))
    return nothing
end


text:
export solvesystem

"""
    (z, τ, κ) = solvesystem(pb::SDCOContext, z::PointPrimalDual{T matT}, τ::T, κ::T) where {T, matT}

    Solves the following system for (dx, dy, ds, dτ, dκ), where A, b, c are held in the pb::SDCOContext structure:

             A(dx)   -b dτ          =  η rP    =  η (bτ - A(x))
    - A*(dy)         +c dτ  -ds     = -η rD    = -η (-A*(y) - s + c τ)
    b⋅dy     -c⋅dx              -dκ =  η rG    =  η (-b⋅y + c⋅x + κ)

    dx + w•ds•w = γμinv(s) - x
    τ dκ + κ dτ = γμ  - τκ

    where μ = (x⋅s + τκ) / (nc + 1)

    **Note** : the first complementarity equation is the symmetrized and pseudo-linearized version of :
        x•ds + s•dx = γμe - x•s
"""
function solvesystem(pb::SDCOContext, x::PointE{T, matT}, y::AbstractArray{T}, s::PointE{T, matT}, τ::T, κ::T, η::T, γ::T) where {T<:Number, matT<:AbstractArray}
    A = pb.A
    c = pb.c
    b = pb.b

    tol = 1e-13

    μ = mu(x, s, τ, κ)
    sinv = inv(s)

    ## Compute right hand sides
    f1 =  η * get_rP(pb, x, τ)
    f2 = -η * get_rD(pb, y, s, τ)
    f3 =  η * get_rG(pb, x, y, κ)
    f4 = product(sinv, γ * μ) - x
    f5 = γ*μ  - τ*κ

    g = NTget_g(pb, x, s)
    w = NTget_w(pb, g)

    M = NTget_M(pb, g)
    factorize(M)

    ## Solving subsystem for dy, dτ
    wcw = hadamard(w, c, w)
    bAwcw = b + evaluate(A, wcw)

    f̂1 = f1 - evaluate(A, f4 + hadamard(w, f2, w))
    f̂2 = f3 + dot(c, f4) + f5 / τ + dot(c, hadamard(w, f2, w))

    sol1 = M \ f̂1
    sol2 = M \ bAwcw

    num   = f̂2 - dot(b, sol1) + dot(c, hadamard(w, evaluate(A, sol1), w))
    denom = dot(b, sol2) - dot(c, hadamard(w, evaluate(A, sol2), w)) + κ/τ + dot(c, wcw)

    dτ = num / denom
    dy = sol1 + sol2*dτ

    @show   norm( evaluate(A, hadamard(w, evaluate(A, dy), w)) - bAwcw * dτ - f̂1 )
    @assert norm( evaluate(A, hadamard(w, evaluate(A, dy), w)) - bAwcw * dτ - f̂1 ) < tol

    @show   norm( dot(b, dy) - dot(c, hadamard(w, evaluate(A, dy), w)) + (κ/τ + dot(c, wcw))*dτ - f̂2 )
    @assert norm( dot(b, dy) - dot(c, hadamard(w, evaluate(A, dy), w)) + (κ/τ + dot(c, wcw))*dτ - f̂2 ) < tol

    ## Deriving full solution

    ds = -evaluate(A, dy) + product(c, dτ) - f2
    dx = f4 - hadamard(w, ds, w)
    dκ = (f5 - κ*dτ) / τ

    # @show norm( evaluate(A, dx) - dτ*b - f1)
    @assert norm( evaluate(A, dx) - dτ*b - f1) < tol
    # @show norm( -evaluate(A, dy) - ds + product(c, dτ) - f2)
    @assert norm( -evaluate(A, dy) - ds + product(c, dτ) - f2) < tol
    # @show norm( dot(b, dy) - dot(c, dx) - dκ - f3)
    @assert norm( dot(b, dy) - dot(c, dx) - dκ - f3) < tol
    # @show norm( dx + hadamard(w, ds, w) - f4)
    @assert norm( dx + hadamard(w, ds, w) - f4) < tol
    # @show norm( τ*dκ + κ*dτ - f5)
    @assert norm( τ*dκ + κ*dτ - f5) < tol

    @assert norm(dot(dx, ds) + dκ*dτ - η*(1-γ-η)*(pb.nc+1)*μ) < tol

    return dx, dy, ds, dτ, dκ
end

function checksystem(pb::SDCOContext{T}, x, y, s, τ, κ, η, γ, dx, dy, ds, dτ, dκ) where {T<:Number}
    μ = mu(x, s, τ, κ)
    sinv = inv(s)
    A = pb.A
    c = pb.c
    b = pb.b

    ## Compute right hand sides
    f1 =  η * get_rP(pb, x, τ)
    f2 = -η * get_rD(pb, y, s, τ)
    f3 =  η * get_rG(pb, x, y, κ)
    f4 = product(ones(c, Dense{T}), γ * μ) - hadamard(x, s)
    f5 = γ*μ  - τ*κ

    tol = 1e-13

    # @show norm( evaluate(A, dx) - dτ*b - f1)
    @assert norm( evaluate(A, dx) - dτ*b - f1) < tol
    # @show norm( -evaluate(A, dy) - ds + product(c, dτ) - f2)
    @assert norm( -evaluate(A, dy) - ds + product(c, dτ) - f2) < tol
    # @show norm( dot(b, dy) - dot(c, dx) - dκ - f3)
    @assert norm( dot(b, dy) - dot(c, dx) - dκ - f3) < tol
    # @show norm( dx + hadamard(w, ds, w) - f4)
    @assert norm( hadamard(dx, s) + hadamard(ds, x) - f4) < tol
    # @show norm( τ*dκ + κ*dτ - f5)
    @assert norm( τ*dκ + κ*dτ - f5) < tol
end

function mu(z::PointPrimalDual{T, matT}, τ::T, κ::T) where {T<:Number, matT<:AbstractArray}
    return mu(x::PointE{T, matT}, s::PointE{T, matT}, τ, κ)
end

function mu(x::PointE{T, matT}, s::PointE{T, matT}, τ::T, κ::T) where {T<:Number, matT<:AbstractArray}
    return (dot(x, s) + τ * κ) / (get_nc(x) + 1)
end

get_rP(pb::SDCOContext, x, τ) = τ * pb.b - evaluate(pb.A, x)
get_rD(pb::SDCOContext, y, s, τ) = -evaluate(pb.A, y) - s + product(pb.c, τ)
get_rG(pb::SDCOContext, x, y, κ) = -dot(pb.b, y) + dot(pb.c, x) + κ

function solve(pb::SDCOContext{T}) where {T<:Number}

    x = ones(pb.c, Dense{T})
    y = zeros(T, pb.m)
    s = ones(pb.c, Dense{T})
    τ = T(1)
    κ = T(1)

    rP0 = get_rP(pb, x, τ)
    rD0 = get_rD(pb, y, s, τ)
    rG0 = get_rG(pb, x, y, κ)

    x0 = ones(pb.c)
    y0 = zeros(T)
    s0 = ones(pb.c)
    τ0 = 1
    κ0 = 1

    μ0 = mu(x, s, τ, κ)

    μ = μ0

    one = ones(pb.c, Dense{T})
    maxit = 5
    it = 0

    μ_hist = T[μ0]
    θ_hist = T[1]
    α_hist = T[]

    ittype = :corr

    print_header(pb)
    print_it(pb, it, x, y, s, τ, κ, rP0, rD0, last(μ_hist), last(θ_hist), 0.01)

    it += 1

    while it < maxit
        printstyled("------ iteration $it, $ittype\n", color=:red)

        μ = mu(x, s, τ, κ)
        println("μ = ", μ)
        println("| x•s - μe | = ", norm(hadamard(x, s) - μ * one))
        println(" | τκ - μ |  = ", norm(τ*κ - μ))

        if ittype == :pred
            η = T(0)
        else
            η = T(1)
        end
        γ = T(1-η)

        ## Solving system
        dx, dy, ds, dτ, dκ = solvesystem(pb, x, y, s, τ, κ, η, γ)
        checksystem(pb, x, y, s, τ, κ, η, γ, dx, dy, ds, dτ, dκ)

        # Property that should always hold
        @assert norm(dot(dx, ds) + dκ*dτ - η*(1-γ-η)*(pb.nc+1)*μ) < 1e-13

        @show norm(dx)
        @show norm(dy)
        @show norm(ds)

        ## Choosing stepsize α
        α = 0.3

        ## Updating variables
        add!(x, α*dx)
        add!(y, α*dy)
        add!(s, α*ds)
        τ += α*dτ
        κ += α*dκ

        ## Variables must remain in their cones
        @assert min(pb, x) > 0
        @assert min(pb, s) > 0
        @assert τ > 0
        @assert κ > 0

        @show min(pb, x), min(pb, s), τ, κ

        rPk = get_rP(pb, x, τ)
        rDk = get_rD(pb, y, s, τ)
        rGk = get_rG(pb, x, y, κ)

        θ = (1-α*η) * last(θ_hist)
        μ_th = (1-α*η) * (1 - α*(1-γ-μ)) * last(μ_hist)
        μ = mu(x, s, τ, κ)

        print_header(pb)
        print_it(pb, it, x, y, s, τ, κ, rPk, rDk, μ, θ, 0.01)

        ## Checking evolution rules hold, logging, it printing

        push!(θ_hist, θ)
        push!(μ_hist, μ)
        push!(α_hist, α)

        # @show norm(μ_th - mu(x, s, τ, κ))

        # @show norm(rPk - θ * rP0)
        # @show norm(rDk - θ * rD0)
        # @show norm(rGk - θ * rG0)

        # @show dot(x0, s) + dot(s0, x) + τ0 * κ + τ * κ0 - (μ/θ + θ*μ0)*(get_nc(x)+1)

        printstyled("Condition 1 (9) : |dx⋅ds + dτ⋅dκ| ≤ O(μ)   :  ", norm(dot(ds, dx) + dτ*dκ), "  ≤  ", μ, "\n", color=:yellow)
        printstyled("Condition 2 (12): Ω(1) ≤ μ/θ ≤ O(1), θ → 0 :  ", θ / μ, "    ", θ, "\n", color=:yellow)

        if ittype == :pred
            ittype = :corr
        else
            ittype = :pred
        end

        it+=1
    end

    @show μ_hist
    @show θ_hist
    @show α_hist

    nothing
end

function print_header(pb::SDCOContext)
    println("it  primerr    dualerr    gaperr       μ          θ          α     it. time")
end

function print_it(pb::SDCOContext, it, x, y, s, τ, κ, rPk, rDk, μ, θ, time)
    primfeaserr = norm(rPk, Inf) / ((1 + norm(pb.b, Inf)) * τ)
    dualfeaserr = norm(rDk, Inf) / ((1 + norm(pb.c, Inf)) * τ)
    gap = abs(dot(pb.c, x) - dot(pb.b, y)) / (τ + max(dot(pb.c, x), dot(pb.b, y)))

    @printf("%2i  %.3e  %.3e  %.3e    %.2e   %.2e   %f\n", it, primfeaserr, dualfeaserr, gap, μ, θ, 1)
end


text:
"""
Finite difference calculation & machine learning of electric potential & electric field from charge
"""

using LinearAlgebra
using Plots
using Random
using Flux
Random.seed!(1)
include("../src/operators.jl")
include("../src/plotutils.jl")
# using Main.EquivariantOperators # end users should omit Main.

# make grid
dims = 2
dx = 0.1
cell = dx * Matrix(I, dims, dims)
rmax = 1.0
grid = Grid(cell, rmax)

# make operators
rmin = 1e-9
ϕ = Op(r -> 1 / (4π * r), rmin, rmax, cell)
E = Op(r -> 1 / (4π * r^2), rmin, rmax, cell; l = 1)
▽ = Op(:▽, cell)

# put dipole charges
ρf = zeros(size(grid))
put!(ρf, grid, [0.5, 0.0],1)
put!(ρf, grid, [-0.5, 0.0],-1)

# calculate fields
Ef = E(ρf)
ϕf = ϕ(ρf)

# test
rvec = [0, 0]
@show get(ϕf, grid, rvec), [0.0]
@show get(Ef, grid, rvec), get(-▽(ϕf), grid, rvec), [-2 / (4π * 0.5^2), 0]

p = []
push!(p, heatmap(ρf',title = "dipole charge"), )
push!(p, heatmap(ϕf', title = "dipole potential"),)
plot(p..., layout = length(p))
vfplot(Ef,grid; title = "dipole electric field")

##
# make neural operators
ϕ_ = Op(Radfunc(), rmin, rmax, cell)
E_ = Op(Radfunc(), rmin, rmax, cell; l = 1)

ps = Flux.params(ϕ_, E_)
function loss()
    remake!(E_)
    remake!(ϕ_)
    global E_f = E_(ρf)
    global ϕ_f = ϕ_(ρf)
    @show l = (nae(E_f, Ef) + nae(ϕ_f, ϕf)) / 2
    l
end

data = [()]
loss()
opt = ADAM(0.1)
Flux.@epochs 10 Flux.train!(loss, ps, data, opt)

## plot
p = []
push!(p, heatmap(ϕ_.kernel, title = "learned potential kernel"))
r = 0:0.01:1
push!(
    p,
    plot(r, ϕ_.radfunc.(r), title = "learned potential kernel radial function"),
)
plot(p..., layout = length(p))

vfplot(E_.kernel,grid; title = "learned E field kernel ")


text:
module TestLeafOptionalArgument

using ComoniconTypes
using ComoniconTargetExpr
using ComoniconTargetExpr: emit_expr, emit_body, emit_norm_body, emit_dash_body
using Test

function foo(a, b=2)
    @test a == 3
end

cmd = Entry(;
    version=v"1.2.0",
    root=LeafCommand(;
        fn=foo,
        name="leaf",
        args=[
            Argument(;name="a", type=Int),
            Argument(;name="b", type=Int, require=false),
        ],
    )
)

eval(emit_expr(cmd))

@testset "test lead optional argument" begin
    @test command_main(["3"]) == 0
    @test command_main(["3", "5"]) == 0
end

end


text:
using .Units
using .PreprocessUnits
using .DependencyParsing
using DependencyParserTest

struct Pipeline
  sentence_splitter::AbstractSentenceSplitter
  tokenizer::AbstractTokenizer
  pos_tagger::AbstractPOSTagger
  dep_parser::AbstractDepParser

  Pipeline(;
    sentence_splitter =  SentenceSplitter(),
    tokenizer = Tokenizer(),
    pos_tagger = POSTagger("tmp\\pos_model"),
    dep_parser = DepParser(DependencyParsing.Settings(), "tmp/success/model_b3000_adamax_fl2_e100_beam_nt_2_best.txt")
  ) = new(sentence_splitter, tokenizer, pos_tagger, dep_parser)
end

function (pipeline::Pipeline)(text)
  completed_trees = Dict();

  mutex = ReentrantLock()

  parse_tree = (sentence_tokens_with_tags, sentence_index) -> begin
    tree = pipeline.dep_parser(sentence_tokens_with_tags)

    lock(mutex) do
      completed_trees[sentence_index] = tree
    end
  end

  sentences = pipeline.sentence_splitter(text)

  Threads.@threads for i in 1:length(sentences)
    tokens_with_tags = pipeline.tokenizer(sentences[i]) |> sentence_tokens -> pipeline.pos_tagger(sentence_tokens)

    parse_tree(tokens_with_tags, i)
  end

  sort(completed_trees |> collect, by = pair -> pair[begin]) |> pairs -> map(pair -> pair[end], pairs)
end

function pipeline_test(connlu_file)
  text = DependencyParserTest.Converting.extract_sentences_from_conllu(connlu_file)
  text = replace(text, r"[éè]" => "e")
  text = replace(text, r"[ô]" => "o")
  pipe = Pipeline()

  trees = pipe(text)

  test_file = "tmp/test_parsed_trees.txt"

  open(test_file, "w") do file
    foreach(trees) do tree
      tree_text = DependencyParser.DependencyParsing.convert_to_string(tree)

      if tree_text != "\n"
        write(file, tree_text)
        write(file, "\n\n")
      end
    end
  end

  conllu_source = DependencyParserTest.Sources.ConlluSource(connlu_file)
  parser_source = DependencyParserTest.Sources.CoreNLPSource(test_file)

  DependencyParserTest.Benchmark.test_accuracy(conllu_source, parser_source)
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
# Automatically generated
struct Params{T}
    z::Vector{Float64}
    footang::T
    g::T
    iff::T
    ihat::T
    irf::T
    ish::T
    ith::T
    k1::T
    k2::T
    k3::T
    k4::T
    k5::T
    k6::T
    k7::T
    k8::T
    lff::T
    lffo::T
    lhat::T
    lmtpxi::T
    lrf::T
    lrff::T
    lrffo::T
    lrfo::T
    lsh::T
    lsho::T
    lth::T
    ltho::T
    ltoexi::T
    mff::T
    mhat::T
    mrf::T
    msh::T
    mth::T
    mtpb::T
    mtpk::T
    rmtpxi::T
    rtoexi::T
    u4::T
    u5::T
    u6::T
    u7::T
    u8::T
    u9::T

end

# initialise with constant values
function Params(footang, g, iff, ihat, irf, ish, ith, k1, k2, k3, k4, k5, k6, k7, k8, lff, lffo, lhat, lmtpxi, lrf, lrff, lrffo, lrfo, lsh, lsho, lth, ltho, ltoexi, mff, mhat, mrf, msh, mth, mtpb, mtpk, rmtpxi, rtoexi)
    z = Vector{Float64}(undef, 827)
    u4 = 0.0
    u5 = 0.0
    u6 = 0.0
    u7 = 0.0
    u8 = 0.0
    u9 = 0.0
    z[380] = g * msh
    z[378] = g * mff
    z[725] = lffo * z[378]
    z[379] = g * mrf
    z[19] = cos(footang)
    z[20] = sin(footang)
    z[68] = lff - lffo
    z[69] = lrf - 0.5lrfo
    z[70] = lsh - lsho
    z[71] = lth - ltho
    z[72] = mhat + 2mff + 2mrf + 2msh + 2mth
    z[74] = (lff * mff + lff * mhat + lffo * mff + 2 * lff * mrf + 2 * lff * msh + 2 * lff * mth) / z[72]
    z[75] = (lrfo * mrf + 2 * lrf * mff + 2 * lrf * mhat + 2 * lrf * mrf + 4 * lrf * msh + 4 * lrf * mth) / z[72]
    z[76] = (lrffo * mrf) / z[72]
    z[77] = (lsh * mff + lsh * mhat + lsh * mrf + lsh * msh + lsho * msh + 2 * lsh * mth) / z[72]
    z[78] = (lth * mff + lth * mhat + lth * mrf + lth * msh + lth * mth + ltho * mth) / z[72]
    z[79] = (mff * z[68]) / z[72]
    z[80] = (lrf * mff + mrf * z[69]) / z[72]
    z[81] = (lsh * mff + lsh * mrf + msh * z[70]) / z[72]
    z[82] = (lth * mff + lth * mrf + lth * msh + mth * z[71]) / z[72]
    z[377] = g * mhat
    z[381] = g * mth
    z[411] = z[74] - lff
    z[412] = 0.5 * z[75] - lrf
    z[413] = z[77] - lsh
    z[414] = z[78] - lth
    z[499] = 0.5 * z[75] - 0.5lrfo
    z[500] = 0.5 * z[76] - 0.5lrffo
    z[524] = 2lrf - z[75]
    z[525] = lff - z[74]
    z[530] = lsh - z[77]
    z[545] = lsh - z[81]
    z[546] = lth - z[82]
    z[547] = lrf - z[80]
    z[548] = lth - z[78]
    z[553] = (lrf - 0.5lrfo) - z[80]
    z[576] = z[19] * z[76]
    z[586] = 2 * z[499] + 2 * z[19] * z[500]
    z[588] = 2 * z[500] + 2 * z[19] * z[499]
    z[596] = 2 * z[553] + 2 * z[19] * z[500]
    z[598] = 2 * z[500] + 2 * z[19] * z[553]
    z[608] = lff * mhat
    z[610] = lrf * mhat
    z[612] = lsh * mhat
    z[614] = lth * mhat
    z[620] = lffo * mff
    z[626] = lff * mrf
    z[628] = lrfo * mrf
    z[630] = lrffo * mrf
    z[672] = mrf * z[69]
    z[683] = lsh * mrf
    z[700] = msh * z[70]
    z[715] = mth * z[71]
    z[727] = z[71] * z[381]
    z[730] = z[70] * z[380]
    z[754] = ihat + 2iff + 2irf + 2ish + 2ith + mff * lffo^2
    z[755] = lrffo^2 + lrfo^2 + 4 * lff^2 + 2 * lrffo * lrfo * z[19]
    z[756] = lff * lrffo
    z[757] = lff * lrfo
    z[758] = lrffo^2 + 4 * z[69]^2
    z[759] = lrffo * z[19] * z[69]
    z[760] = lff^2 + lrf^2 + lsh^2 + lth^2
    z[761] = lff * lsh
    z[762] = lff * lth
    z[763] = lrf * lth
    z[764] = lff * lrf
    z[765] = lrf * lsh
    z[766] = lsh * lth
    z[768] = lrffo * z[19]
    z[769] = lrfo * z[19]
    z[770] = lrffo * z[20]
    z[771] = lrfo * z[20]
    z[772] = z[19] * z[69]
    z[773] = z[20] * z[69]
    z[775] = iff + irf + ish + ith
    z[776] = lrffo^2
    z[777] = z[69]^2
    z[782] = iff + irf + ish + ith + mff * lffo^2
    z[787] = iff + irf + ish
    z[788] = lsh * z[69]
    z[789] = lrffo * lsh
    z[794] = iff + irf + ish + mff * lffo^2
    z[795] = lff^2 + lrf^2 + lsh^2
    z[800] = iff + irf
    z[801] = (lrffo^2 + 4 * z[69]^2) - 4 * lrffo * z[19] * z[69]
    z[806] = iff + irf + mff * lffo^2
    z[807] = lff^2 + lrf^2

    return Params(z, footang, g, iff, ihat, irf, ish, ith, k1, k2, k3, k4, k5, k6, k7, k8, lff, lffo, lhat, lmtpxi, lrf, lrff, lrffo, lrfo, lsh, lsho, lth, ltho, ltoexi, mff, mhat, mrf, msh, mth, mtpb, mtpk, rmtpxi, rtoexi, u4, u5, u6, u7, u8, u9)
end


text:
#__precompile__()

#using ArgCheck
#using Knet
#using Pkg; Pkg.build("Knet")
# Remove these later

#=
Julia implementation of Mask R-CNN
=#

#include("rpn.jl")
#include("fpn.jl")
#inclide("resnet.jl")

struct Mask_RCNN
    rpn
    fpn
    anchors
end

function config(config = nothing)
    # Copying most of the config from matterport/mask_RCNN
    # Define resnet
    # Define FPN

    # Can change name depending on experiement
    NAME = nothing

    # Number of GPUs to use. When using only a CPU, this needs to be set to 1.
    GPU_COUNT = 1

    # Number of images to train with on each GPU. A 12GB GPU can typically handle 2 images of 1024x1024
    # REMOVE ME BEFORE PUBLISHING - currently have IMAGES_PER_GPU set to 1 because my GPU (1070ti) has 8GB memory
    IMAGES_PER_GPU = 1

    # Number of training steps per epoch
    STEPS_PER_EPOCH = 1000

    # Number of validation steps to run at the end of every training epoch
    VALIDATION_STEPS = 50

    # Backbone network architecture
    # This supports resnet50, resnet101, and resnet152
    BACKBONE = "resnet101"

    # The strides of each layer of the FPN. These values are based on the default resnet101 backbone
    BACKBONE_STRIDES = [4, 8, 16, 32, 64]

    # The shape of the backbone
    BACKBONE_SHAPES = [256 256; 128 128; 64 64; 32 32; 166 16;]

    # Size of the fully-connected layers in the classifier
    FPN_CLASSIF_FC_LAYERS_SIZE = 1024

    # Size of the top-down layers used to build the feature pyramid
    TOP_DOWN_PYRAMID_SIZE = 256

    # Number of classification classes (including background) -> overwite this in implementation
    NUM_CLASSES = 1

    # Length of square anchor side in pixels
    RPN_ANCHORS_SCALES = [32, 64, 128, 256, 512]

    # Ratios of anchors at each cell (width/height)
    # A value of 1 represents a square anchor, and 0.5 is a wide anchor
    RPN_ANCHOR_RATIOS = [0.5, 1, 2]

    # Anchor stride
    # If 1, then anchors are created for each cell in the backbone feature map
    # If 2, then anchors are created for every other cell, and so on.
    RPN_ANCHOR_STRIDE = 1

    # Non-max suppression threshold tot filter RPN proposals
    # You can increase this during training to generate more proposals
    RPN_NMS_THRESHOLD = 0.7

    # How many anchors per image to use for RPN training
    RPN_TRAIN_ANCHORS_PER_IMAGE = 256

    # Define pool size
    POOL_SIZE = 7
end


text:
combustion = ce"C2H4 + O2 = CO2 + 2 H2O"
redox = ce"H{+} + Cr2O7{2-} + C2H5OH = Cr{3+} + CO2 + H2O"
cerational = ChemEquation{Rational}("1//2 H2 + O2 → H2O")
cefloat = ChemEquation{Float64}("0.33 N3 + 0.5 O2 = NO")

@testset "equationmatrix" begin
    @test equationmatrix(combustion) == [
        2  0  1  0
        4  0  0  2
        0  2  2  1
    ]
    @test equationmatrix(redox) == [
        1   0  6  0  0  2
        0   2  0  1  0  0
        0   7  1  0  2  1
        0   0  2  0  1  0
        1  -2  0  3  0  0
    ]
    @test equationmatrix(cerational) == [
        2//1  0//1  2//1
        0//1  2//1  1//1
    ]
    @test equationmatrix(cefloat) == [
        3.0  0.0  1.0
        0.0  2.0  1.0
    ]
    @test equationmatrix(ce"151 H2 + 32 O2 = 7 H2O") ==
        equationmatrix(ce"5 H2 + O2 = 5 H2O")
end

@testset "balancematrix" begin
    # [:,:] required to reshape Vector to Nx1 Array
    @test balancematrix(combustion) == [1//1; 3//1; -2//1; -2//1][:,:]
    @test balancematrix(redox, fractions=true) == [16//11; 2//11; 1//11; -4//11; -2//11; -1//1][:,:]
    @test balancematrix(cerational) == [1//1; 1//2; -1//1][:,:]
    @test balancematrix(cefloat) == [0.3333333333333333; 0.5; -1.0][:,:]
    @test balancematrix(ce"H2 + O2 = H2O") == balancematrix(ce"2 H2 + O2 = 2 H2O")
end

@testset "balance" begin
    # Original equation should be unchanged
    balanced_combustion = balance(combustion)
    @test combustion == combustion

    # Examples from https://www.webqc.org/balance.php
    equations = [
        "Fe + Cl2 = FeCl3" =>
            "2 Fe + 3 Cl2 = 2 FeCl3",
        "KMnO4 + HCl = KCl + MnCl2 + H2O + Cl2" =>
            "2 KMnO4 + 16 HCl = 2 KCl + 2 MnCl2 + 8 H2O + 5 Cl2",
        "K4Fe(CN)6 + H2SO4 + H2O = K2SO4 + FeSO4 + (NH4)2SO4 + CO" =>
            "K4Fe(CN)6 + 6 H2SO4 + 6 H2O = 2 K2SO4 + FeSO4 + 3 (NH4)2SO4 + 6 CO",
        "C6H5COOH + O2 = CO2 + H2O" =>
            "2 C6H5COOH + 15 O2 = 14 CO2 + 6 H2O",
        "K4Fe(CN)6 + KMnO4 + H2SO4 = KHSO4 + Fe2(SO4)3 + MnSO4 + HNO3 + CO2 + H2O" =>
            "10 K4Fe(CN)6 + 122 KMnO4 + 299 H2SO4 = 162 KHSO4 + 5 Fe2(SO4)3 + 122 MnSO4 + 60 HNO3 + 60 CO2 + 188 H2O",
        "Cr2O7{-2} + H{+} + {-} = Cr{+3} + H2O" =>
            "Cr2O7{-2} + 14 H{+} + 6 e = 2 Cr{+3} + 7 H2O",
        "S{-2} + I2 = I{-} + S" =>
            "S{-2} + I2 = 2 I{-} + S",
        "PhCH3 + KMnO4 + H2SO4 = PhCOOH + K2SO4 + MnSO4 + H2O" =>
            "5 PhCH3 + 6 KMnO4 + 9 H2SO4 = 5 PhCO2H + 3 K2SO4 + 6 MnSO4 + 14 H2O",
        "CuSO4*5H2O = CuSO4 + H2O" =>
            "CuSO4*5H2O = CuSO4 + 5 H2O",
        "Zn(s) + O2(g) = ZnO(s)" => "2 Zn + O2 = 2 ZnO"
    ]
    for equation ∈ equations
        @test balance(ChemEquation(equation[1])) == ChemEquation(equation[2])
    end

    @test balance(cerational) == ChemEquation{Rational}("H2 + 1//2 O2 = H2O")
    @test balance(cefloat) == ChemEquation{Float64}("0.3333333333333333 N3 + 0.5 O2 = NO")
    @test balance(combustion, fractions=true) == ChemEquation{Rational}("1//2 C2H4 + 3//2 O2 = CO2 + H2O")

    @test_throws ErrorException balance(ce"H2 + O = H + O")
    @test_throws ErrorException balance(ce"H2 + CO = H2O")
end


text:
@doc raw"""
    energy_balance(canonical::Canonical,
                        initial_conditions::Vector{InitialCondition},
                        efficiency_data::Tuple{Vector{String}, Vector{InOut}},
                        cons_name::Symbol,
                        var_names::Tuple{Symbol, Symbol, Symbol})

Constructs multi-timestep constraint from initial condition, efficiency data, and variable tuple

# Constraints

If t = 1:

``` varenergy[name, 1] == initial_conditions[ix].value + varin[name, 1]*eff_in*fraction_of_hour - varout[name, 1]*fraction_of_hour/eff_out ```

If t > 1:

``` varenergy[name, t] == varenergy[name, t-1] + varin[name, t]*eff_in*fraction_of_hour - varout[name, t]*fraction_of_hour/eff_out ```

# LaTeX

`` x^{energy}_1 == x^{energy}_{init} + frhr \eta^{in} x^{in}_1 - \frac{frhr}{\eta^{out}} x^{out}_1, \text{ for } t = 1 ``

`` x^{energy}_t == x^{energy}_{t-1} + frhr \eta^{in} x^{in}_t - \frac{frhr}{\eta^{out}} x^{out}_t, \forall t \geq 2 ``

# Arguments
* canonical::Canonical : the canonical model built in PowerSimulations
* initial_conditions::Vector{InitialCondition} : for time zero 'varenergy'
* efficiency_data::Tuple{Vector{String}, Vector{InOut}} :: charging/discharging efficiencies
* cons_name::Symbol : name of the constraint
* var_names::Tuple{Symbol, Symbol, Symbol} : the names of the variables
- : var_names[1] : varin
- : var_names[2] : varout
- : var_names[3] : varenergy

"""
function energy_balance(canonical::Canonical,
                        initial_conditions::Vector{InitialCondition},
                        efficiency_data::Tuple{Vector{String}, Vector{InOut}},
                        cons_name::Symbol,
                        var_names::Tuple{Symbol, Symbol, Symbol})

    time_steps = model_time_steps(canonical)
    resolution = model_resolution(canonical)
    fraction_of_hour = Dates.value(Dates.Minute(resolution))/60
    name_index = efficiency_data[1]

    varin = get_variable(canonical, var_names[1])
    varout = get_variable(canonical, var_names[2])
    varenergy = get_variable(canonical, var_names[3])

    constraint = add_cons_container!(canonical, cons_name, name_index, time_steps)

    for (ix, name) in enumerate(name_index)
        eff_in = efficiency_data[2][ix].in
        eff_out = efficiency_data[2][ix].out

        constraint[name, 1] = JuMP.@constraint(canonical.JuMPmodel,
                                   varenergy[name, 1] == initial_conditions[ix].value + varin[name, 1]*eff_in*fraction_of_hour
                                                    - (varout[name, 1])*fraction_of_hour/eff_out)

    end

    for t in time_steps[2:end], (ix, name) in enumerate(name_index)
        eff_in = efficiency_data[2][ix].in
        eff_out = efficiency_data[2][ix].out

        constraint[name, t] = JuMP.@constraint(canonical.JuMPmodel,
                                   varenergy[name, t] == varenergy[name, t-1] + varin[name, t]*eff_in*fraction_of_hour
                                                    - (varout[name, t])*fraction_of_hour/eff_out)
    end

    return

end


text:
macro stridedloops(N, itersym, dimsym, args...)
    _stridedloops(N, itersym, dimsym, args...)
end
function _stridedloops(N::Int, itersym::Symbol, dimsym::Symbol, args...)
    mod(length(args),3)==1 || error("Wrong number of arguments")
    body = args[end]
    ex = Expr(:escape, body)
    for i=1:3:length(args)-1
        ex=Cartesian.sreplace!(ex,args[i],Cartesian.inlineanonymous(args[i], 1))
    end
    for dim = 1:N
        itervar = Cartesian.inlineanonymous(itersym, dim)
        dimvar = Cartesian.inlineanonymous(dimsym, dim)
        preargs = Any[]
        postargs = Any[]
        for i=1:3:length(args)-1
            indnew = Cartesian.inlineanonymous(args[i], dim)
            start = (dim < N ? Cartesian.inlineanonymous(args[i], dim+1) : args[i+1])
            step = Cartesian.inlineanonymous(args[i+2], dim)
            push!(preargs,:($(esc(indnew)) = $(esc(start))))
            push!(postargs,:($(esc(indnew)) += $(esc(step))))
        end
        preexpr=Expr(:block,preargs...)
        postexpr=Expr(:block,postargs...)
        ex = quote
            $preexpr
            for $(esc(itervar)) = 1:$(esc(dimvar))
                $ex
                $postexpr
            end
        end
    end
    if N == 0 # even a zero-dimensional array has one element and thus requires a single run
        preargs = Any[]
        for i=1:3:length(args)-1
            indnew = Cartesian.inlineanonymous(args[i], 1)
            start = args[i+1]
            push!(preargs,:($(esc(indnew)) = $(esc(start))))
        end
        preexpr=Expr(:block,preargs...)
        ex = quote
            $preexpr
            $ex
        end
    end
    ex
end

macro gentracekernel(N1,N2,order,alpha,A,beta,C,startA,startC,odims,cdims,ostridesA,cstridesA,ostridesC)
    _gentracekernel(N1,N2,order,alpha,A,beta,C,startA,startC,odims,cdims,ostridesA,cstridesA,ostridesC)
end
function _gentracekernel(N1::Int,N2::Int,order::Symbol,alpha::Symbol,A::Symbol,beta::Symbol,C::Symbol,
    startA::Symbol,startC::Symbol,odims::Symbol,cdims::Symbol,ostridesA::Symbol,cstridesA::Symbol,ostridesC::Symbol)
    ex=quote
        local indA1, indA2, indC
        local i,j
        if $(esc(order))==0
            local gamma
            gamma=$(esc(beta))
            @stridedloops($N1, i, $(esc(cdims)), indA1, $(esc(startA)), $(esc(cstridesA)), begin
                @stridedloops($N2, j, $(esc(odims)), indA2, indA1, $(esc(ostridesA)), indC, $(esc(startC)), $(esc(ostridesC)), begin
                    @inbounds $(esc(C))[indC]=gamma*$(esc(C))[indC]+$(esc(alpha))*$(esc(A))[indA2]
                end)
                gamma=one($(esc(beta)))
            end)
        else
            @stridedloops($N2, j, $(esc(odims)), indA1, $(esc(startA)), $(esc(ostridesA)), indC, $(esc(startC)), $(esc(ostridesC)), begin
                local localC
                @inbounds localC=$(esc(beta))*$(esc(C))[indC]
                @stridedloops($N1, i, $(esc(cdims)), indA2, indA1, $(esc(cstridesA)), @inbounds localC+=$(esc(alpha))*$(esc(A))[indA2])
                @inbounds $(esc(C))[indC]=localC
            end)
        end
    end
    ex
end

macro gencontractkernel(N1,N2,N3,order,alpha,A,conjA,B,conjB,beta,C,startA,startB,startC,odimsA,odimsB,cdims,ostridesA,cstridesA,ostridesB,cstridesB,ostridesCA,ostridesCB)
    _gencontractkernel(N1,N2,N3,order,alpha,A,conjA,B,conjB,beta,C,startA,startB,startC,odimsA,odimsB,cdims,ostridesA,cstridesA,ostridesB,cstridesB,ostridesCA,ostridesCB)
end
function _gencontractkernel(N1::Int,N2::Int,N3::Int,order::Symbol,
    alpha::Symbol,A::Symbol,conjA::Symbol,B::Symbol,conjB::Symbol,beta::Symbol,C::Symbol,
    startA::Symbol,startB::Symbol,startC::Symbol,odimsA::Symbol,odimsB::Symbol,cdims::Symbol,
    ostridesA::Symbol,cstridesA::Symbol,ostridesB::Symbol,cstridesB::Symbol,ostridesCA::Symbol,ostridesCB::Symbol)
    ex=quote
        local indA1, indA2, indB1, indB2, indC1, indC2, gamma
        # we still have to implement other orders
        if $(esc(order))==3 # k,j,i
            @stridedloops($N3, k, $(esc(cdims)), indA1, $(esc(startA)), $(esc(cstridesA)), indB1, $(esc(startB)), $(esc(cstridesB)), begin
                gamma=$(esc(beta))
                @stridedloops($N2, j, $(esc(odimsB)), indB2, indB1, $(esc(ostridesB)), indC1, $(esc(startC)), $(esc(ostridesCB)), begin
                    @inbounds localB=($(esc(conjB))=='C' ? conj($(esc(B))[indB2]) : $(esc(B))[indB2])
                    @stridedloops($N1, i, $(esc(odimsA)), indA2, indA1, $(esc(ostridesA)), indC2, indC1, $(esc(ostridesCA)), begin
                        @inbounds localA=($(esc(conjA))=='C' ? conj($(esc(A))[indA2]) : $(esc(A))[indA2])
                        @inbounds localC=gamma*$(esc(C))[indC2]
                        @inbounds $(esc(C))[indC2]=localC+$(esc(alpha))*localA*localB
                    end)
                end)
                gamma=one($(esc(beta)))
            end)
        else # j,i,k
            @stridedloops($N2, j, $(esc(odimsB)), indB1, $(esc(startB)), $(esc(ostridesB)), indC1, $(esc(startC)), $(esc(ostridesCB)), begin
                @stridedloops($N1, i, $(esc(odimsA)), indA1, $(esc(startA)), $(esc(ostridesA)), indC2, indC1, $(esc(ostridesCA)), begin
                    @inbounds localC=$(esc(beta))*$(esc(C))[indC2]
                    @stridedloops($N3, k, $(esc(cdims)), indA2, indA1, $(esc(cstridesA)), indB2, indB1, $(esc(cstridesB)), begin
                        @inbounds localA=($(esc(conjA))=='C' ? conj($(esc(A))[indA2]) : $(esc(A))[indA2])
                        @inbounds localB=($(esc(conjB))=='C' ? conj($(esc(B))[indB2]) : $(esc(B))[indB2])
                        localC+=$(esc(alpha))*localA*localB
                    end)
                    @inbounds $(esc(C))[indC2]=localC
                end)
            end)
        end
    end
    ex
end


text:
module ImgTools

    import Images
    import FileIO
    import FixedPointNumbers: N0f8
    import Colors: RGB
    import Plots
    import Plots: AbstractPlot
    import DataFileNames: dfname

    include("sgif.jl")
    include("grid.jl")
    include("ticks.jl")
    include("sfig.jl")
    include("lfig.jl")

    export sfig, lfig, sgif

end


text:
#
#  Hénon-Heiles
#

function HenonHam(u,params)

    @inbounds begin
    q = view(u,1:2)
    p = view(u,3:4)

    T=zero(eltype(u))
    V=zero(eltype(u))

    T=1/2*(p[1]^2+p[2]^2)
    V=1/2*(q[1]^2+q[2]^2)+q[1]^2*q[2]-1/3*q[2]^3

    return(T+V)
    end

end

#
#  Dynamical Pro
#

function dotq(dq,q,p,params,t)
    @inbounds begin
    dq[1] = p[1]
    dq[2] = p[2]
    end
end

function dotp(dp,q,p,params,t)
    @inbounds begin
    dp[1] = -q[1]-2*q[1]*q[2]
    dp[2] = -q[2]-q[1]^2+q[2]^2
    end
end

#
#  Second Order Problem
#

function f2nd!(ddu,du,u,p,t)

  @inbounds begin
  q = view(u,1:2)   # x

  ddu[1] = -q[1]-2*q[1]*q[2]
  ddu[2] = -q[2]-q[1]^2+q[2]^2

  end
end


text:
module SpaceMethods

export AbstractSpaceMethod, AbstractDGMethod, odefun!

"""
    AbstractSpaceMethod

Supertype for spatial discretizations
"""
abstract type AbstractSpaceMethod end

"""
    AbstractDGMethod <: AbstractSpaceMethod

Supertype for discontinuous Galerkin spatial discretizations
"""
abstract type AbstractDGMethod <: AbstractSpaceMethod end

"""
    odefun!(disc::AbstractSpaceMethod, dQ, Q, t)

Evaluates the right-hand side of the spatial discretization defined by `disc` at
time `t` with state `Q`. The result is added into `dQ`. Namely, the
semi-discretization is of the form
``
    Q̇ = F(Q, t)
``
and after the call `dQ += F(Q, t)`

!!! note

    There is no generic implementation of this function. This must be
    implemented for each subtype of `AbstractSpaceMethod`
"""
odefun!(m::AbstractSpaceMethod, dQ, Q, t) =
throw(MethodError(odefun!, (m, dQ, Q, t)))

end

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
println("loading packages... ")

using ArgParse
using OrthographicProjection

println("packages OK")

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
		"source"
			help = "A text file with Potree directories list or a single Potree directory"
			required = true
		"--output", "-o"
            help = "Output file"
			required = true
		"--bbin"
            help = "Bounding box as 'x_min y_min z_min x_max y_max z_max' or Potree JSON volume model"
			required = true
		"--po"
            help = "Projection plane: XY+, XZ+, YZ+"
			arg_type = String
            default = "XY+"
		"--quote"
			help = "Distance of plane from origin"
			arg_type = Float64
		"--thickness"
			help = "Section thickness"
			arg_type = Float64
		"--ucs"
			help = "Path to UCS JSON file. If not provided is the identity matrix"
			arg_type = String
		"--epsg"
			help = "EPSG code"
			arg_type = Int
    end

    return parse_args(s)
end

function main()
    args = parse_commandline()

	OrthographicProjection.flushprintln("== params ==")
    for (arg,val) in args
        OrthographicProjection.flushprintln("$arg  =>  $val")
    end

	bbin = args["bbin"]
	output = args["output"]
	PO = args["po"]
	txtpotreedirs = args["source"]
	q = args["quote"]
	thickness = args["thickness"]
	ucs = args["ucs"]
	epsg = args["epsg"]

	b = tryparse.(Float64,split(bbin, " "))
	if length(b) == 6
		#bbin = (hcat([b[1],b[2],b[3]]),hcat([b[4],b[5],b[6]]))
		bbin = Common.AABB(b[4],b[1],b[5],b[2],b[6],b[3])
	end

	if isnothing(ucs)
		ucs = Matrix{Float64}(Common.I,4,4)
	else
		ucs = FileManager.ucs2matrix(ucs)
	end

	coordsystemmatrix = OrthographicProjection.PO2matrix(PO,ucs)
	model = Common.getmodel(bbin)
	aabb = Common.boundingbox(model[1])

	if !isnothing(q) && !isnothing(thickness)
		origin = Common.inv(ucs)[1:3,4]
		model = Common.getmodel(Common.inv(coordsystemmatrix), q, thickness, aabb; new_origin = origin)
	end

	seg = OrthographicProjection.segment(txtpotreedirs, output, model; epsg = epsg)
	proj_folder = splitdir(output)[1]
	FileManager.successful(seg, proj_folder)
end

@time main()


text:
using BenchmarkTools
using ACE

function dosomething!(A, x)
   @inbounds begin
      A[1] = 1; A[2] = x
      for n = 3:length(A)
         A[n] = 2 * x * A[n-1] - A[n-2]
      end
   end
   return A
end

pool = ACE.

A1 = zeros(30)

po@btime dosomething!($A, 0.123)
@btime dosomething!($A1, 0.123)
@btime dosomething!($A2, 0.123)


text:
module YaoLang

using LinearAlgebra
using YaoAPI

include("runtime/locations.jl")
include("runtime/generic_circuit.jl")

module Compiler

using TimerOutputs
const to = TimerOutput()
timings() = (TimerOutputs.print_timer(to); println())
enable_timings() = (TimerOutputs.enable_debug_timings(Compiler); return)

using ExprTools
using IRTools
using IRTools.Inner
using YaoAPI

using IRTools.Inner:
    IR, Block, BasicBlock, Statement, Variable, block, blocks, tab, print_stmt, printargs
using Base.Meta: ParseError

import ..YaoLang
using YaoLang: AbstractLocations, merge_locations, Locations, CtrlLocations, Circuit, PrimitiveCircuit

include("compiler/parse.jl")
include("compiler/ir.jl")
include("compiler/print.jl")
include("compiler/codegen.jl")
include("compiler/compiler.jl")
include("compiler/reflection.jl")
include("compiler/primitive.jl")
include("compiler/utils.jl")
include("compiler/kwdefs.jl")
include("compiler/validation.jl")
include("compiler/trace.jl")

function __init__()
    TimerOutputs.reset_timer!(to)
    # not sure why this doesn't work inside the module
    IRTools.Inner.printers[:quantum] = function (io, ex)
        get(printers, ex.args[1], print)(io, ex)
    end
end
end

using .Compiler
export @device, @primitive, @code_yao, @quantum
include("runtime/primitives.jl")

end # module


text:
module OndaEDF

using Base.Iterators
using Dates, UUIDs
using Onda, EDF
using StatsBase
using TimeSpans
using Tables
using Tables: rowmerge

# can be dropped if we drop Onda<0.14
sample_type(x) = isdefined(Onda, :sample_type) ? Onda.sample_type(x) : x.sample_type

include("standards.jl")

include("import_edf.jl")
export edf_to_onda_samples, edf_to_onda_annotations, store_edf_as_onda
export extract_channels_by_label, edf_signals_to_samplesinfo

include("export_edf.jl")
export onda_to_edf

include("prettyprint_diagnostic_info.jl")

end # module


text:
using Distributions

module RandomProcesses
	using Distributions

	export ChineseRestaurantProcess, CRP
	export ChineseRestaurantTableProcess, CRT

	include("crp.jl")
	include("crt.jl")
end


text:
addprocs()
using MAT
include("include.jl")

TopN = 10

# Read Data
DATA = matread("yahoo.mat")
TrainSet = DATA["TrainSet"]
Holdout = DATA["Holdout"]
UW = DATA["SampledUnwatched"]
n, m = size(TrainSet)

# Read Item Model
W = readItemModel("example.model",m)

# Build RecWalk Model
P = RecWalk(TrainSet,W,0.005)

NDCG = SharedArray{Float64}(n)
RR = SharedArray{Float64}(n)
HR = SharedArray{Float64}(n)

#Base Item Model
PI = TrainSet*W
@sync @parallel for user = 1:n
	HR[user], RR[user], NDCG[user] = Single_HR_RR_NDCG(PI[user,:], vcat(Holdout[user], UW[:,user]), TopN)
end
println("Base Item Model:  HR = $(mean(HR))  ARHR=$(mean(RR))  NDCG=$(mean(NDCG))")

# RecWalk - K-Step
K = 7
@sync @parallel for user = 1:n
    ru  = sparse(reshape(P[user,:], 1, m+n))
    [ru  *= P for step=2:K]
    HR[user], RR[user], NDCG[user] = Single_HR_RR_NDCG(ru[n+1:end], vcat(Holdout[user], UW[:,user]), TopN)
end
println("RecWalk K-Step:   HR = $(mean(HR))  ARHR=$(mean(RR))  NDCG=$(mean(NDCG))")

# RecWalk - PR
eta = 0.7
PI = inv(full(I-eta*P)) # due to the small size of the example data the recwalk ppr vectors can be computed in batch.
PI = PI[1:n,n+1:end]
@sync @parallel for user = 1:n
	HR[user], RR[user], NDCG[user] = Single_HR_RR_NDCG(PI[user,:], vcat(Holdout[user], UW[:,user]), TopN)
end
println("RecWalk PR:       HR = $(mean(HR))  ARHR=$(mean(RR))  NDCG=$(mean(NDCG))")

max_stars_repo_name: nikolakopoulos/RecWalk

max_issues_repo_name: nikolakopoulos/RecWalk

max_forks_repo_name: nikolakopoulos/RecWalk


text:
# Add parent::Int here to avoid overwriting same method in JuMP. Should be safe.
function JuMP.parseNLExpr_runtime(m::JuMP.Model, x::JuMP.Variable, tape, parent::Int, values)
    #JuMP.__last_model[1] = x.m  #__last_model is not defined in JuMP
    if x.m === m
        push!(tape, JuMP.NodeData(ReverseDiffSparse.VARIABLE, x.col, parent))
    else
        othermap = getStructure(m).othermap
        if haskey(othermap, x)
            newx = othermap[x]
        else
            newx = JuMP.Variable(m,JuMP.getlowerbound(x),JuMP.getupperbound(x),:Cont, JuMP.getname(x))
            othermap[x] = newx
        end
        push!(tape, JuMP.NodeData(ReverseDiffSparse.VARIABLE, newx.col, parent))
    end
    nothing
end


text:
using Soss
using MeasureTheory

m = @model n begin
    α ~ Normal()
    β ~ Normal()
    x ~ Normal() |> iid(n)
    σ ~ Exponential(λ=1)
    y ~ For(x) do xj
        Normal(α + β * xj, σ)
    end
    return y
end

rand(m(3))
rand(m(3),5)
simulate(m(3), 5)
mysim = simulate(m(3), 1000)
mysim[1]

mytrace = mysim.trace

using TupleVectors
@with mytrace begin
        ŷ = α .+ β .* x
        r = y - ŷ
        (;ŷ, r)
end


text:
"A module to plot Ulam spirals"
module UlamSpiral

import PyPlot
using Primes
import Base: position

export position, text_spiral, point_spiral, plot_quadratic

"""
    position(n[, start_n=1])

Compute (x,y) location in the lattice spiral of the integer n, assuming the spiral
starts at the integer start_n (which takes default value 1).

# Arguments
- `n::Integer`: the integer to find the position of, N.B. n >= start_n
- `start_n::Integer=1`: the starting integer of the spiral

# Examples
```jldoctest
julia> UlamSpiral.position(23, 1)
(0, -2)
```
"""
function position( n::Integer, start_n::Integer=1 )
    if n<start_n
        throw(DomainError())
    elseif n==start_n
        return (0, 0)
    end

    j = n-start_n
    k = Int(ceil( floor(sqrt( j )) / 2 )) # how far out in the spiral
    m = Int( j - (2*(k-1)+1)^2 ) # how far along the arm
    length = Int( 2*k )
    quarter = Int( floor(m / length) )
    m2 = m - length*quarter
    if quarter==0
        d_x = k
        d_y =  Int( m2 - length/2 + 1 )
    elseif quarter==1
        d_x = -Int( m2 - length/2 + 1 )
        d_y = k
    elseif quarter==2
        d_x = -k
        d_y = -Int( m2 - length/2 + 1 )
    elseif quarter==3
        d_x =  Int( m2 - length/2 + 1 )
        d_y = -k
    end
    # println(i, ", ", k, ", ", m, ", ", length, ", ", quarter, ", ", d_x, ", ", d_y)
    return (d_x, d_y)
end

"""
    text_spiral( [start_n] )

Write numbers in an Ulam spiral starting from start_n.

# Arguments
- `start_n::Integer=1`: the starting integer of the spiral

# Keyword Arguments
- `fontsize=20`:        the size of the font to write everything
- `color_prime="red"`:  the color to use for prime numbers
- `color_comp="blue"`:  the color to use for composite numbers (note 1 is compositite by convention)
- `width=10`:           the width of the square used (which determines how many numbers to plot)
- `fig_n=1`:            the figure on which to plot everything
- `clear="on"`:         whether to clear the figure before starting

# Examples
```jldoctest
julia> text_spiral()
PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x1221f3dd0>)

julia> text_spiral( 100; fontsize=16)
PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x1221f3dd0>)

julia> text_spiral(fig_n=2, color_comp="white", color_prime="black", width=20, fontsize=10)
PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x122136e10>)

```
"""
function text_spiral(start_n::Integer=1;
                     fontsize=20, color_prime="red", color_comp="blue", width=10, fig_n=1, clear="on")
    c_x = Int(floor(width/2)) # centre of spiral
    c_y = Int(floor(width/2))
    stop_n = start_n + width^2 - 1;

    f = figure(fig_n)
    if clear=="on"
        clf()
    end
    hold(true)
    for n=start_n : stop_n
        (d_x, d_y) = UlamSpiral.position(n, start_n)
        if isprime(n)
            PyPlot.text(c_x + d_x, c_y + d_y, @sprintf("%d", n), fontsize=fontsize, color=color_prime)
        else
            PyPlot.text(c_x + d_x, c_y + d_y, @sprintf("%d", n), fontsize=fontsize, color=color_comp)
         end
    end
    PyPlot.xlim(0, width+1)
    PyPlot.ylim(0, width+1)
    return f
end

"""
    point_spiral( [start_n] )

Plot points at numbers in an Ulam spiral starting from start_n.

# Arguments
- `start_n::Integer=1`: the starting integer of the spiral

# Keyword Arguments
- `marker="o":     the marker to use
- `markersize=7`:  the size of the markers being plotted
- `color="red"`:   the color to use (this routine only plots prime numbers)
- `width=10`:      the width of the square used (which determines how many numbers to plot)
- `fig_n=1`:       the figure on which to plot everything
- `clear="on"`:    whether to clear the previous figure

# Examples
```jldoctest
julia> point_spiral()
PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x1221f3dd0>)

julia> point_spiral( 1000; markersize=10)
PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x1221f3dd0>)

julia> point_spiral(fig_n=2, markersize=1, color="black", width=200)
PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x122136e10>)

julia> text_spiral(fig_n=2, color_comp="white", color_prime="black", width=20, fontsize=10)
julia> point_spiral(fig_n=2, markersize=3, marker="s", color="green", width=20, clear="off")
PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x122136e10>)

```
"""
function point_spiral(start_n::Integer=1;
                      marker="o", markersize=7, color="red", width=100, fig_n=1, clear="on")
    c_x = Int(floor(width/2)) # centre of spiral
    c_y = Int(floor(width/2))
    stop_n = start_n + width^2 - 1;

    f = figure(fig_n)
    if clear=="on"
        clf()
    end
    hold(true)
    for n=primes(start_n, stop_n)
        (d_x, d_y) = UlamSpiral.position(n, start_n)
        PyPlot.plot(c_x + d_x, c_y + d_y; marker=marker, markersize=markersize, color=color)
    end
    PyPlot.xlim(0, width+1)
    PyPlot.ylim(0, width+1)
    return f
end

"""
    plot_quadratic(a, b, c)

Highlight points on the quadratic ax^2 + bx + c. Presumes you have already
run point_spiral with the same width parameter.

# Arguments
- `a::Integer = 4`:  quadratic coefficient
- `b::Integer = -2`: linear coefficient
- `c::Integer = 41`: constant coefficient

# Keyword Arguments
- `start_n=1`           the starting integer of the spiral
- `marker="o":          the marker to use
- `markersize=7`:       the size of the markers being plotted
- `color_prime="blue"`: the color to use for prime numbers
- `color_comp="white"`: the color to use for compositive numbers
- `width=10`:              the width of the square used (should be the same as previous plot)
- `fig_n=1`:              the figure on which to plot everything (should be the same as previous plot)

# Examples
```jldoctest
julia> point_spiral(fig_n=2, markersize=1, color="black", width=200);
julia> plot_quadratic(; fig_n=2, width=200);
julia> plot_quadratic(4, 2, 41; fig_n=2, width=200, color_prime="green");
```
"""
function plot_quadratic(a::Integer=4, b::Integer=-2, c::Integer=41;
                        start_n=1, width=10, fig_n=1,
                        marker="o", markersize=7, color_prime="red", color_comp="white")
    c_x = Int(floor(width/2)) # centre of spiral
    c_y = Int(floor(width/2))
    stop_n = start_n + width^2 - 1;
    rad = 0
    x = start_n - 1
    while rad <= width/2
        x = nextprime(x+1)
        n = a*x^2 + b*x + c
        if n>0
            (d_x, d_y) = UlamSpiral.position(n, start_n);
            rad = max( abs(d_x), abs(d_y) )
            plot(c_x + d_x, c_y + d_y, marker=marker, color=color_prime, markersize=markersize)
        end
    end
    # at present not plotting composite numbers at all
end

end # module


text:
using Escher
using FactCheck
using Compat

include("macros.jl")
include("interop.jl")
FactCheck.exitstatus()


text:
# the loop version of loop
function compact!(x)
    @inbounds for i in 1:3
        if x[i] == 0
            for j in i+1:4
                if x[j] != 0
                    tmp = x[j:4]
                    x[j:4] .= 0
                    x[i:i+4-j] .= tmp
                    break;
                end
            end
        end
    end
    x
end

function move_loop!(x)
    # firstly, compact x so that there are no gaps
    # however cells are not yet merged
    compact!(x)

    # now the data is compact; we combine it
    @inbounds if (x[1] != 0) && (x[1] == x[2])
        x[1] += 1
        if (x[3] != 0)  && (x[3] == x[4])
            x[2], x[3], x[4] = x[3] + 1, 0 , 0
        else
            x[2], x[3], x[4] = x[3], x[4], 0
        end
    elseif (x[2] != 0) && (x[2] == x[3])
        x[2] += 1
        x[3], x[4] = x[4], 0
    elseif (x[3] != 0) && (x[3] == x[4])
        x[3] += 1
        x[4] = 0
    end

    x
end


text:
# Define a (custom) Standard Normal distribution, for illustrative puspose.
struct StandardNormal <: Distributions.ContinuousUnivariateDistribution end
Distributions.logpdf(::StandardNormal, x::Real) = -(x ^ 2 + log(2 * pi)) / 2
Distributions.rand(rng::AbstractRNG, ::StandardNormal) = randn(Random.GLOBAL_RNG)


text:
include("Facto.jl")
include("testPrimitive.jl")
include("testFacto.jl")

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
include("miscellaneous.jl")
include("def.jl")
include("util_iter.jl")
include("util_macros.jl")
include("error.jl")
include("runtime.jl")
include("agent.jl")
include("signal.jl")
include("aql.jl")
include("memory.jl")
include("queue.jl")
include("code.jl")
include("config.jl")
include("extensions.jl")


text:
module DrugRepositioning

import Base:show,size,sum,+,eltype
using LinearAlgebra
using DelimitedFiles
using Random:shuffle!
using BSON
using CUDA
using CUDA:CuArray
using Statistics
using Distances:pairwise,CosineDist,PreMetric,Metric,SemiMetric
using Zygote
using Flux
using Flux:gpu
import Flux.Data:_getobs,_nobs
using Plots
using DataFrames
import Plots:plot

DisType=Union{PreMetric,SemiMetric,Metric}

include("dataset/preprocess.jl")
include("dataset/dataset.jl")
include("experiment/crossvalidatation.jl")
include("experiment/metrics.jl")
include("experiment/compare.jl")
include("model/io.jl")
include("model/model.jl")
include("model/train.jl")

end # module

max_stars_repo_name: yeruoforever/DrugRepositioning

max_issues_repo_name: yeruoforever/DrugRepositioning

max_forks_repo_name: yeruoforever/DrugRepositioning


text:
# NOTE: This does what its supposed to do but I don't think its
# a good idea to use it actually.
# NOTE: I just copied this from wikipedia, not sure how it actually works
function hilbert_ordering(d::Int, n::Int)

    t = d
    x = 0.0
    y = 0.0

    for i=1:n
        s = 2^(i-1)
        rx = 1 & div(t, 2)
        ry = 1 & xor(t, rx)

        if (ry == 0)
            if (rx == 1)
                x = s - 1 - x
                y = s - 1 - y
            end

            r = Int[y, x]
            x = r[1]
            y = r[2]
        end

        x += s * rx
        y += s * ry
        t = div(t, 4)
    end

    return Int[x, y]
end


text:
using QuadGK
using SpecialFunctions
using Roots
using RCall

function Gammadist(gamma)
	return ((adap.be^adap.al)/SpecialFunctions.gamma(adap.al))*(gamma^(adap.al-1))*exp(-adap.be*gamma)
end

function pip0(gamma)
	U = 4*adap.theta_f*adap.Lf/(2.0*adap.NN)
	R = 2*adap.Lf*adap.rho/(2.0*adap.NN)
	return Gammadist(gamma)*exp(-(Gammadist(gamma)*U/(2.0*adap.NN))/(gamma/(adap.NN+0.0)+R/(2.0*adap.NN)))
end

function intpip0()

	f(gam) = pip0(gam)
	return QuadGK.quadgk(f,0,1000)[1]
end

function calcBGam(L,alpha,beta,theta)
	# not sure whats going on here, seems to overestimate or under sometimes
	u = 2*theta/(2.0*adap.NN)
	r = adap.rho/(2.0*adap.NN)

	a = -1.0*(2^(alpha+1))*exp(L*adap.NN*r*beta)*L*adap.N*((1.0/(L*adap.N*r))^(1-alpha))*u*(beta^alpha)
	b =  convert(Float64,SpecialFunctions.gamma_inc(1-alpha,L*adap.NN*r*beta,1)[2])
	#fudge = 1-gamInt.cdf(0.2,a=adap.al2,scale=1./adap.be2)
	#fudge = 1.
	fudge = 0.25

	c = exp(a*b*fudge)

	return c
end

function set_theta_f_gam()
	i(theta) = calcBGam(adap.Lf,adap.al2,adap.be2,theta)-adap.B
	# theta_f  = fsolve(lambda theta ,0.0000001)
	theta_f  = Roots.find_zero(i,0.0000001)
	adap.theta_f = theta_f
end

function calcB(L,theta)
	t = -1.0*adap.gam_neg/(adap.NN+0.)
	u = 2.0*theta/(2.0*adap.NN)
	r = adap.rho/(2.0*adap.NN)

	#return u*t/((t+r*L)^2)
	#return u*t/(2*(t+(1.-np.exp(-2*r*L))/2.)^2)

	# Nordborg models
	#####return u/(2*r*t*(1+((1.-np.exp(-2*r*L))/2.)*(1-t)/t)^2)
	#return u/(t*(1+((1-np.exp(-2.0*r*L))/2.)*(1-t)/t)^2)
	#return u/(t*(1+r*L*(1-t)/t)^2)
	#####return u/(t*(1+(np.exp(-2*r*L)/2)/t)^2)
end

# Br() needs L and theta, not working right now
function get_B_vals()
	ret = Array{Float64}(undef,30)
	for i in 20:50
		L = convert(Int64,round(1.3^i,digits=0))
		ret[i] = (Br(L),L)
	end
	return ret
end

"""
	Estimating and plotting MAP using locfit and ggplot2 in R. It assume your folder contains the posterior estimated through ABCreg
"""
function plotMap(;analysis::String,output::String)

	try
		@eval using RCall
		@eval R"""library(ggplot2);library(abc)"""

		out = filter(x -> occursin("post",x), readdir(analysis,join=true))
		out = filter(x -> !occursin(".1.",x),out)

		open(x) = Array(CSV.read(GZip.open(x),DataFrame,header=false))
		posteriors = open.(out)

		maxp = DataFrame(Array{Float64,2}(undef,size(posteriors,1),5),[:aw,:as,:a,:gamNeg,:shape])
		R"""getmap <- function(df){
				temp = as.data.frame(df)
			    d <-locfit(~temp[,1],temp);
			    map<-temp[,1][which.max(predict(d,newdata=temp))]
			}"""
		getmap(x) = rcopy(R"""matrix(apply($x,2,getmap),nrow=1)""")
		tmp = getmap.(posteriors)
		maxp = DataFrame(vcat(tmp...),[:aw,:as,:a,:gamNeg,:shape])

		al  = maxp[:,1:3]
		gam  = maxp[:,4:end]
		p = rcopy(R"""al = $al
			dal = melt(al)
			pal = ggplot(dal) + geom_density(aes(x=value,fill=variable),alpha=0.5) + scale_fill_manual(values=c('#30504f', '#e2bd9a', '#ab2710'))
			ggsave(pal,filename=paste0($output))
			""")
		return(maxp)
	catch
		println("Please install R, ggplot2 and abc in your system before execute this function")
	end
end


text:
module TestUniform

using StatsBase
using Test
using TreeParzen
import TreeParzen: API, Trials

N = 10_000
COUNTMAX = 1300
COUNTMIN = 700

samples = [TreeParzen.Resolve.node(HP.Uniform(:u, 0.0, 10.0), TreeParzen.Trials.ValsDict()) for i in 1:N]

@test minimum(samples) > 0
@test maximum(samples) < 10
h = fit(StatsBase.Histogram, samples, nbins=10).weights
@test all(COUNTMIN .< h)
@test all(h .< COUNTMAX)

end
true


text:
using StructuresKit

#Test shear and moment calculations for a cantilever beam.

#Ix Iy Ixy J Cw
sectionProperties = [(100.0, 100.0, 0.0, 10.0, 10.0)]

#E  ν
materialProperties = [(29500, 0.30)]

#ax ay
loadLocation = [(0.0, 0.0)]

#kx kϕ
springStiffness = [(0.0, 0.0)]

#member information
#L dL SectionProperties MaterialProperties LoadLocation SpringStiffness
memberDefinitions = [(30*12,3,1,1,1,1)]

#end boundary conditions
#type=1 u''=v''=ϕ''=0 (simply supported), type=2 u'=v'=ϕ'=0  (fixed), type=3 u''=v''=ϕ''=u'''=v'''=ϕ'''=0 (free end, e.g., a cantilever)
endBoundaryConditions = [2 3]

#supports
#location where u=v=ϕ=0
supports = [0.0]

#load  qx   qy
uniformLoad = (0.0, 5.0/1000)

z, u, v, ϕ, beamProperties = Beam.solve(memberDefinitions, sectionProperties, materialProperties, loadLocation, springStiffness, endBoundaryConditions, supports, uniformLoad)

Mxx = InternalForces.moment(z, beamProperties.dm, -v, beamProperties.E, beamProperties.Ix)
Vyy = InternalForces.shear(z, beamProperties.dm, -v, beamProperties.E, beamProperties.Ix)

maxMxx = minimum(Mxx)
maxMxx_AISC = -uniformLoad[2]*z[end]^2/2

maxVyy = maximum(Vyy)
maxVyy_AISC = uniformLoad[2]*z[end]

MxxError= abs((maxMxx - maxMxx_AISC))/ maxMxx_AISC
VyyError= abs((maxVyy - maxVyy_AISC))/ maxVyy_AISC

#accept 1% error from numerical solution
@test MxxError <= 0.01
@test VyyError <= 0.01


text:
using FilePathsBase
using FilePathsBase: /, join
using Base.Filesystem
using Dates: datetime2unix
using JLSO
using Mmap
using Test

using FilePathsBase.TestPaths

include("testpkg.jl")

@testset "FilePathsBase" begin
    include("mode.jl")
    include("buffer.jl")
    include("system.jl")

    @static if Sys.isunix()
        # Test that our weird registered path works
        ps = PathSet(TestPkg.posix2test(tmpdir()) / "pathset_root"; symlink=true)

        @testset "$(typeof(ps.root))" begin
            testsets = [
                test_registration,
                test_show,
                test_cmd,
                test_parse,
                test_convert,
                test_components,
                test_indexing,
                test_iteration,
                test_parents,
                test_descendants_and_ascendants,
                test_join,
                test_splitext,
                test_basename,
                test_filename,
                test_extensions,
                test_isempty,
                test_normalize,
                test_canonicalize,
                test_relative,
                test_absolute,
                test_isdir,
                test_isfile,
                test_stat,
                test_filesize,
                test_modified,
                test_created,
                test_cd,
                test_readpath,
                test_walkpath,
                test_read,
                test_write,
                test_mkdir,
                test_cp,
                test_mv,
                test_symlink,
                test_touch,
                test_tmpname,
                test_tmpdir,
                test_mktmp,
                test_mktmpdir,
                test_download,
                test_issocket,
                # These will also all work for our custom path type,
                # but many implementations won't support them.
                test_isfifo,
                test_ischardev,
                test_isblockdev,
                test_ismount,
                test_isexecutable,
                test_isreadable,
                test_iswritable,
                test_chown,
                test_chmod,
            ]

            # Run all of the automated tests
            test(ps, testsets)
        end
    end
end


text:
using JuLIP, ForwardDiff, StaticArrays, BenchmarkTools

# simple EAM potential
f(R) = sqrt( 1.0 + sum( exp(-norm(r)) for r in R ) )

nneigs = 30
R0 = [ @SVector rand(3) for n = 1:nneigs ]

f_fd(R) =  ForwardDiff.gradient( T -> f(vecs(T)),  mat(R)[:] ) |> vecs

function f_d(R::AbstractVector{<:JVec})
   ∇f = zeros(JVecF, length(R))
   ρ̄ = sum( exp(-norm(r)) for r in R )
   dF = 0.5 * (1.0 + ρ̄)^(-0.5)
   return [ dF * (- exp(-norm(r))) * r/norm(r)  for r in R ]
end

@assert f_d(R0) ≈ f_fd(R0)

@btime f($R0);
@btime f_d($R0);
@btime f_fd($R0);

# f1(R::AbstractMatrix) =
#       sqrt( 1.0 + sum( exp.(-sqrt.(sum(abs2, R, 1))) ) )
# f1_rd(R::Matrix) =
#       reshape(
#             ReverseDiff.gradient( S -> f1(reshape(S, 3, length(S) ÷ 3)), R[:] ),
#             size(R) )
# R1 = mat(R0)
# f1_rd(R1)
#
# @btime f1($R1);
# @btime f1_rd($R1);

##

module M
   using ForwardDiff
   myf(x) = x^2+x+1
   myg(x) = π   # deliberately incorrect gradient for testing

   function myf(d::ForwardDiff.Dual{T}) where T
      x = ForwardDiff.value(d)
      @show d
      # y = ( myf(x), myg(x) )
      ForwardDiff.Dual{T}( myf(x), myg(x) * ForwardDiff.partials(d) )
   end
end

using ForwardDiff
M.myf(1.0)
ForwardDiff.derivative(M.myf, 1.0)

##

module S
   using ForwardDiff, Dierckx

   modelspline() = Spline1D([0.0, 0.5, 1.0], rand(3); bc = "zero", k = 2)
   ev(s::Spline1D, x) = s(x)
   ev_d(s::Spline1D, x) = Dierckx.derivative(s, x, 1)
   ev_dd(s::Spline1D, x) = Dierckx.derivative(s, x, 2)

   function ev(s::Spline1D, d::ForwardDiff.Dual{T}) where T
      x = ForwardDiff.value(d)
      ForwardDiff.Dual{T}( s(x), ev_d(s, x) * ForwardDiff.partials(d) )
   end

   function ev_d(s::Spline1D, d::ForwardDiff.Dual{T}) where T
      x = ForwardDiff.value(d)
      ForwardDiff.Dual{T}( ev_d(s, x), ev_dd(s, x) * ForwardDiff.partials(d) )
   end
   # Dierckx.derivative(V.spl, r, nu)

   function compf(s, X)
      [ sum( ev_d(s, x)^p for x in X )  for p = 0:3 ]
   end

   function J_compf(s, X)
      [ (p * ev_d(s, X[i])^(p-1) * ev_dd(s, X[i])) for p=0:3, i=1:length(X) ]
   end
end

using ForwardDiff
s = S.modelspline()
ForwardDiff.derivative( x -> S.ev(s, x), 0.5) == S.ev_d(s, 0.5)
ForwardDiff.derivative( x -> S.ev_d(s, x), 0.5) == S.ev_dd(s, 0.5)

X = rand(5)
J = ForwardDiff.jacobian(X_ -> S.compf(s, X_), X)
J1 = S.J_compf(s, X)
J1 ≈ J

##

using JuLIP, Test, Printf, LinearAlgebra, SparseArrays, ForwardDiff
using JuLIP.Potentials, JuLIP.Testing

data = joinpath(dirname(pathof(JuLIP)), "..", "data") * "/"
eam_Fe = JuLIP.Potentials.EAM(data * "pfe.plt", data * "ffe.plt", data * "F_fe.plt")

##

at = set_pbc!( bulk(:Fe, cubic = true), false ) * 2
rattle!(at, 0.1)
nlist = neighbourlist(at, cutoff(eam_Fe))
_, Rs = neigs(nlist, 1)
Xs = [ [zero(JVecF)]; Rs ]
x0 = mat(Xs)[:]

function _conf2env(x)
   Xs = vecs(x)
   return [ Xs[j] - Xs[1] for j = 2:length(Xs) ]
end

_dV2conf(dV) = mat([[- sum(dV)]; dV])[:]

fsite = X -> JuLIP.evaluate(eam_Fe, _conf2env(X))
fsite_d = X -> _dV2conf( JuLIP.evaluate_d(eam_Fe, _conf2env(X)) )

@test fsite(x0) == JuLIP.evaluate(eam_Fe, Rs)
@test JuLIP.Testing.fdtest(fsite, fsite_d, x0)

# ------

fsite_dd = X -> ForwardDiff.jacobian(fsite_d, X)
JuLIP.Testing.fdtest_hessian(fsite_d, fsite_dd, x0)

H = fsite_dd(x0)
Hp = zeros(JMatF, 11, 11)
for i = 1:11, j = 1:11
   Hp[i, j] = JMatF( H[ i*3 .+ (1:3), j*3 .+ (1:3) ] ... )
end

H1 = JuLIP.Potentials.evaluate_dd(eam_Fe, Rs)
H1 ≈ Hp

##

at = set_pbc!( bulk(:Fe, cubic = true), false ) * 2
rattle!(at, 0.1)
set_calculator!(at, eam_Fe)
# hessian(eam_Fe, at)
println(@test fdtest(eam_Fe, at, verbose=true))
# println(@test fdtest_hessian( x->gradient(at, x), x->hessian(at, x), dofs(at) ))
# fdtest_hessian(x->gradient(at, x), x->hessian(at, x), dofs(at), verbose=true
#       )

H = hessian(eam_Fe, at)
fdtest_hessian(x->gradient(at, x), x->hessian(at, x), dofs(at), verbose=true
      )

x = dofs(at)
U = 0.1 * rand(JVecF, length(at))
u = mat(U)[:]
g1 = H * u
h = 1e-3
g2 = (gradient(at, x+1e-3 * u) - gradient(at, x-1e-3*u)) / (2*h)
g1 - g2

##


text:
function reset_distribution(p::LaserTagPOMDP, b::ParticleCollection, a, o)
    # warn("Resetting Particle Filter Distribution")
    rob = first(particles(b)).robot
    nextrob = LaserTag.add_if_inside(p.floor, rob, LaserTag.ACTION_DIRS[a])
    if o == LaserTag.C_SAME_LOC
        return ParticleCollection{LaserTag.LTState}([LaserTag.LTState(nextrob, nextrob, false)])
    else
        return LaserTag.LTInitialBelief(nextrob, p.floor)
    end
end

max_possible_weight(pomdp::LaserTagPOMDP, a::Int, o) = 0.0

new_particle(pomdp::LaserTagPOMDP, a::Int, o) = error("tried to generate a new particle (shouldn't get here)")

# max_possible_weight(pomdp::LaserTagPOMDP, a::Int, o::Float64) = max(1.0, pdf(Normal(0.0, pomdp.return_std), 0.0))


text:
module MLFlowLogger

using PyCall
using UUIDs
using Requires
using FileIO
using ColorTypes

export MLFLogger

const mlflow = PyNULL()

function __init__()
    copy!(mlflow, pyimport("mlflow"))
    @require ImageIO="82e4d734-157c-48bb-816b-45c225c6df19" begin
        include("image.jl")
        @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
            include("plots.jl")
        end
    end
end

using Base.CoreLogging: CoreLogging, AbstractLogger, LogLevel, Info, handle_message, shouldlog, min_enabled_level, catch_exceptions

mutable struct MLFLogger <: AbstractLogger
    client::PyObject
    run::PyObject
    step_increment::Int
    global_step::Int
    min_level::LogLevel
end

function MLFLogger(; min_level=Info, step_increment=1, start_step=0, experiment_name=nothing, tracking_uri=nothing, kwargs...)
    if isnothing(tracking_uri)
        client = mlflow.tracking.MlflowClient()
    else
        client = mlflow.tracking.MlflowClient(tracking_uri)
    end

    expid = nothing
    if haskey(ENV, "MLFLOW_EXPERIMENT_ID")
        expid = ENV["MLFLOW_EXPERIMENT_ID"]
        exp = client.get_experiment(expid)
        experiment_name = exp.name
    end

    if isnothing(experiment_name)
        experiment_name = string(UUIDs.uuid4())
    end
    if isnothing(expid)
        expid = client.create_experiment(experiment_name)
    end

    if haskey(ENV, "MLFLOW_RUN_ID")
        runid = ENV["MLFLOW_RUN_ID"]
        run = client.get_run(runid)
    else
        run = client.create_run(expid)
    end

    MLFLogger(client, run, step_increment, start_step, min_level)
end

increment_step!(lg::MLFLogger, Δ_Step) = lg.global_step += Δ_Step

add_tag!(lg::MLFLogger, tag::String) = lg.client.set_tag(lg.run.info.run_id, "tag", tag)
add_tag!(lg::MLFLogger, key::String, value::String) = lg.client.set_tag(lg.run.info.run_id, key, value)

"""
    function log_metric(lg::CLogger, name::AbstractString, value::Real; step::Int=nothing, epoch::Int=nothing)

Logs general scalar metrics.
"""
function log_metric(lg::MLFLogger, key::AbstractString, value; timestamp=nothing, step=nothing)
    lg.client.log_metric(lg.run.info.run_id, key, value, timestamp=timestamp, step=step)
end

function log_param(lg::MLFLogger, key::AbstractString, value)
    lg.client.log_param(lg.run.info.run_id, key, value)
end

"""
    function log_artifact(lg::MLFLogger, local_path, artifact_path)

Log a local file or directory as an artifact of the currently active run.

- `local_path`: Path to the file to write
- `artifact_path`: If provided, the directory to write to.
"""
function log_artifact(lg::MLFLogger, local_path, artifact_path=nothing)
    lg.client.log_artifact(lg.run.info.run_id, local_path, artifact_path)
end

"""
    function log_image(lg::MLFLogger, key::AbstractString, obj::AbstractString; step=nothing)

Log a local image file as an artifact of the currently active run. If current `step` provided, adds
it to the artifact directory path.
"""
function log_image(lg::MLFLogger, key::AbstractString, obj::AbstractString; step=nothing)
    if !ispath(obj)
        @warn "$obj not a path to an image"
        return
    end
    if !isnothing(step)
        key = key * "/$step"
    end
    key = key * "/" * basename(obj)
    log_artifact(lg, obj, key)
end

const log_parameter = log_param

CoreLogging.catch_exceptions(lg::MLFLogger) = false

CoreLogging.min_enabled_level(lg::MLFLogger) = lg.min_level

CoreLogging.shouldlog(lg::MLFLogger, level, _module, group, id) = true

function preprocess(name, val::T, data) where {T}
    if isstructtype(T)
        fn = logable_propertynames(val)
        for f=fn
            prop = getfield(val, f)
            preprocess(name*"/$f", prop, data)
        end
    else
        # If we do not know how to serialize a type, then
        # it will be simply logged as text
        push!(data, name=>val)
    end
    data
end

"""
    logable_propertynames(val::Any)
Returns a tuple with the name of the fields of the structure `val` that
should be logged to Comet.ml. This function should be overridden when
you want Comet.ml to ignore some fields in a structure when logging
it. The default behaviour is to return the  same result as `propertynames`.
See also: [`Base.propertynames`](@ref)
"""
logable_propertynames(val::Any) = propertynames(val)

## Default unpacking of key-value dictionaries
function preprocess(name, dict::AbstractDict, data)
    for (key, val) in dict
        # convert any key into a string, via interpolating it
        preprocess("$name/$key", val, data)
    end
    return data
end

# Split complex numbers into real/complex pairs
preprocess(name, val::Complex, data) = push!(data, name*"/re"=>real(val), name*"/im"=>imag(val))

# Handle standard float metrics
process(lg::MLFLogger, name::AbstractString, obj::Real, step::Int) = log_metric(lg, name, obj; step=step)

function CoreLogging.handle_message(lg::MLFLogger, level, message, _module, group, id, file, line; kwargs...)
    i_step = lg.step_increment # :log_step_increment default value

    if !isempty(kwargs)
        data = Vector{Pair{String,Any}}()
        # ∀ (k-v) pairs, decompose values into objects that can be serialized
        for (key,val) in pairs(kwargs)
            # special key describing step increment
            if key == :log_step_increment
                i_step = val
                continue
            end
            preprocess(strip(message*"/$key", ['(', ')', '{', '}', '!']), val, data)
        end
        iter = increment_step!(lg, i_step)
        for (name, val) in data
            process(lg, name, val, iter)
        end
    end
end

Base.show(io::IO, lg::MLFLogger) = begin
	str  = "MLFLogger(\"run=$(lg.run.info.run_id)\", min_level=$(lg.min_level), "*
		   "current_step=$(lg.global_step))"
    Base.print(io, str)
end

Base.close(lg::MLFLogger) = lg.client.set_terminated(lg.run.info.run_id)

end