text:
using NOAu111
using Test
using JuLIP
using PyCall
using ASE
using NQCBase
using NQCModels
using Unitful, UnitfulAtomic

build = pyimport("ase.build")

slab = build.fcc111("Au", size=(2,4,4), vacuum=10.0)
no = build.molecule("NO")
build.add_adsorbate(slab, no, 5, "ontop")

at = JuLIP.Atoms(ASEAtoms(slab))
rattle!(at, 0.5)

@testset "H11" begin
    @test JuLIP.Testing.fdtest(NOAu111.AuNCoupling(), at, verbose=false)
    @test JuLIP.Testing.fdtest(NOAu111.H11DiabaticElement(), at, verbose=false)
    @test JuLIP.Testing.fdtest(NOAu111.ImagePotential(at), at, verbose=false)
end

@testset "Components" begin
    @test JuLIP.Testing.fdtest(NOAu111.AuAu(), at, verbose=false)
    @test JuLIP.Testing.fdtest(NOAu111.H00(), at, verbose=false)
    @test JuLIP.Testing.fdtest(NOAu111.H01(), at, verbose=false)
    @test JuLIP.Testing.fdtest(NOAu111.H11(at), at, verbose=false)
end

@testset "Final model" begin
    R = austrip.(Matrix(hcat(at.X...))u"Å")
    c = NQCBase.PeriodicCell(austrip.(cell(at)' .* u"Å"))

    model = NOAu(chemical_symbols(at), c, R)

    @test JuLIP.Testing.fdtest(x->potential(model, x), x->derivative(model, x), R, verbose=false)
end


text:
###############################################################################
##### Weierstrass elliptic functions
###############################################################################

"""
Set of Weierstrass elliptic functions parameters for computational issues
"""
struct Weierstrass <: EllipticFunction
    ω1            ::  ComplexF64
    ω3            ::  ComplexF64
    e1            ::  ComplexF64
    g2            ::  ComplexF64
    g3            ::  ComplexF64
    η1            ::  ComplexF64
    # (\\frac{π θ_3(0) θ_4(0)}{2ω_1})^2
    # coefficient in ℘-function evaluation
    ℘_factor      ::  ComplexF64
    σ_factor      ::  ComplexF64
    # θ-function data
    θ             ::  Theta
end

"""
  WeierstrassData(; ω1::ComplexF64, ω3::ComplexF64)

Construct WeierstrassData object from lattice generators ω1 and ω3
"""
function Weierstrass(;
        ω1::ComplexF64,
        ω3::ComplexF64,
        θ ::Theta
        )
    #####
    tolerance = eps(imag(ω3))

    th2 = theta(θ, th_k=2, d_n=0, z=0.0im)
    th4 = theta(θ, th_k=4, d_n=0, z=0.0im)
    th24 = th2^4
    th44 = th4^4
    e_factor = π^2 / (12 * ω1^2)
    e1 = e_factor * (th24 + 2th44)
    e2 = e_factor * (th24 - th44)
    e3 = -e_factor * (2th24 + th44)

    g2 = -4 * (e2*e3 + e3*e1 + e1*e2)
    g3 = 4 * e1 * e2 * e3

    d3th1 = theta(θ, th_k=1, d_n=3, z=0.0im)
    d1th1 = theta(θ, th_k=1, d_n=1, z=0.0im)

    η1 = -π^2 / (12ω1) * d3th1 / d1th1
    η3 = (η1 * ω3 - 1im*π/2) / ω1

    th3 = theta(θ, th_k=3, d_n=0, z=0.0im)
    ℘_factor = (π * th3 * th4 / (2ω1))^2

    σ_factor = 2ω1 / (π*d1th1)

    return Weierstrass(ω1, ω3, e1, g2, g3, η1, ℘_factor, σ_factor, θ)
end

"""
weierstrass_normalized!(; z::ComplexF64, norm_factor::Float64, w::WeierstrassData,
        upper_term::Int64, output::ComplexOffsetVector)

Compute sequence of normalized weierstrass elliptic functions derivatives

# Return
output[-2] = ``\\wp^{(-2)}(z)``,
output[-1] = ``\\wp^{(-1)}(z)``,
output[n]  = ``\\frac{r^{n+2}}{(n+1)!} \\wp^{(n)}(z)``, n ≥ 0.
"""
function weierstrass_normalized!(;
        z           ::ComplexF64,
        norm_factor ::Float64,
        w           ::Weierstrass,
        upper_term  ::Int64,
        output      ::ComplexOffsetVector
        )
    #####
    ω1 = w.ω1
    e1 = w.e1
    η1 = w.η1
    g2 = w.g2
    ℘_factor = w.℘_factor
    σ_factor = w.σ_factor

    u = (π / (2ω1)) * z
    th1, th2, dth1, dth2 = theta(w.θ, u)
    th1_norm = th1 / norm_factor

    ζ = (π/(2ω1)) * (dth1/th1_norm)
    ζ += norm_factor * (η1 / ω1) * z

    ℘ = ℘_factor * (th2 / th1_norm)^2
    ℘ += norm_factor^2 * e1

    ℘_prime = th1 * th2 * dth2 - dth1 * th2 * th2
    ℘_prime /= th1_norm^3
    ℘_prime *= (π / ω1) * ℘_factor

    ℘_pp = 6*℘*℘ - 0.5 * g2 * norm_factor^4

    output[-1] = -ζ
    output[0] = ℘             # /0!
    output[1] = ℘_prime       # /1!
    output[2] = ℘_pp / 2.0    # /2!

    # wp_{n+2}(z) = frac{6}{(n+1)(n+2)} sum_{k=0}^n wp_k(z) wp_{n-k}(z),
    # n >= 1.
    for n in 1 : upper_term-2
        output[n+2] = 0.0im
        for k in 0 : n
            output[n+2] += output[k] * output[n-k]
        end
        output[n+2] *= 6.0 / ((n+1) * (n+2))
    end

    # for n in 1 : upper_term
    #     output[n] /= n+1
    # end

    output[-2] = -(η1 * (z^2) / (2ω1) + log(σ_factor * th1))

end

function expand_at_pole!(;
        output           ::ComplexOffsetVector,
        weierstrass_data ::Weierstrass,
        upper_term       ::Int64
        )
    #####
    if upper_term < 2
        error("Too few terms of expansion")
    end

    g2 = weierstrass_data.g2
    g3 = weierstrass_data.g3

    n_c = div(upper_term+2, 2)
    c = OffsetVector{ComplexF64}(undef, 2:n_c)

    c[2] = g2 / 20.0

    if n_c >= 3
        c[3] = g3 / 28.0
    end

    for n in 4 : n_c
        c[n] = 0
        for m in 2 : n-2
            c[n] += c[m] * c[n-m]
        end
        c[n] *= 3 / ((2n+1)*(n-3))
    end

    for n in 0 : upper_term
        output[n] = 0
    end
    for n in 2 : n_c
        output[2n-2] = c[n]
    end

    output[-1] = 0
    output[-2] = 0
end

function term_expansion_normalized!(;
        output         ::ComplexOffsetVector,
        dict_of_series ::RComplex2OffsetMatrix,
        delta          ::RationalComplex,
        n_derivative   ::Int64,
        upper_term     ::Int64,
        r_source       ::Float64,
        r_dest         ::Float64
        )

    # При решении задачи нужно вычислять разложения функций Вейерштрасса
    # в противоположных точках. Чтобы не делать двойную работу, значения вычисляются
    # только для z с положительной действительной частью, а для противоположного -
    # нечетные члены рядов берутся с противоположным знаком.
    if real(delta) < 0
        r = -r_dest  # (-1)^n r^n = (-r)^n
        delta = -delta
    else
        r = r_dest
        # delta = delta
    end

    # Массив разложений в ряд Тейлора производных ℘-функции Вейерштрасса
    # в точке delta.
    α = dict_of_series[delta]

    if delta != 0
        factor = (r_source / abs(delta))^(n_derivative+2)
    else
        factor = r_dest^(n_derivative+2)
    end

    for n in 0 : upper_term
        output[n] = factor * α[n, n_derivative]
        factor *= r
    end
end

ext: jl


text:
@safetestset "multivariate normal" begin
    using DifferentialEvolutionMCMC, Test, Random, Distributions
    using MCMCChains, LinearAlgebra
    Random.seed!(50514)

    # number of variables
    n_μ = 30
    # number of observations per variable
    n_d = 100
    # μ parameters
    μs = fill(0.0, n_μ)
    # data
    data = rand(MvNormal(μs, 1.0 * I), n_d)

    # function for initial values
    function sample_prior()
        μ = rand(Normal(0, 1), n_μ)
        σ = rand(truncated(Cauchy(0, 1), 0, Inf))
        return as_union([μ,σ])
    end

    # returns prior log likelihood
    function prior_loglike(μ, σ)
        LL = 0.0
        LL += sum(logpdf.(Normal(0, 1), μ))
        LL += logpdf(truncated(Cauchy(0, 1), 0, Inf), σ)
        return LL
    end

    # likelihood function
    function loglike(data, μs, σ)
        return sum(logpdf(MvNormal(μs, σ^2 * I), data))
    end

    # upper and lower bounds of parameters
    bounds = ((-Inf,Inf),(0.0,Inf))
    # parameter names
    names = (:μ,:σ)

    # model object
    model = DEModel(;
        sample_prior,
        prior_loglike,
        loglike,
        data,
        names
    )

    # DEMCMC sampler
    de = DE(;
        sample_prior,
        bounds,
        sample = resample,
        burnin = 5000,
        n_initial = (n_μ + 1) * 4,
        Np = 3,
        n_groups = 1,
        θsnooker = 0.1,
    )
    # sample from the posterior distribution
    n_iter = 50_000
    chains = sample(model, de, MCMCThreads(), n_iter, progress=true)
    sds = describe(chains)[1][1:n_μ,:std]
    means = describe(chains)[1][1:n_μ,:mean]
    @test all(x -> isapprox(x, 0.1; atol = .01), sds)
    @test all(x -> isapprox(x, 0.0; atol = .2), means)
    @test std(means) ≈ 0.1 atol = .01
end


text:
module RPlot

    using RCall

    include("graphics.jl")
    include("ggplot.jl")

end # module


text:
####### TO MOVE
import ProximalOperators: prox, prox!, gradient, gradient!

using RecursiveArrayTools
function prox(f::ProximableFunction, x::ArrayPartition, args...)
  y, fy = prox(f, x.x, args...)
  return ArrayPartition(y), fy
end

prox!(y::ArrayPartition, f::ProximableFunction, x::ArrayPartition, args...) =
prox!(y.x, f, x.x, args...)

function gradient(f::ProximableFunction, x::ArrayPartition)
  gradx, fy = gradient(f, x.x)
  return ArrayPartition(gradx), fy
end

gradient!(y::ArrayPartition, f::ProximableFunction, x::ArrayPartition) =
gradient!(y.x, f, x.x)
###############


text:
using JuMP
using CPLEX
include("instances/exemple-lecture-graphe.jl")

function Benders(file::String, b::Int64)
  n, m, adj, demande = readGraph(file::String)
  #### probleme maitre ####
  # Create the model
  m0 = Model(CPLEX.Optimizer)

  ## Variables
  @variable(m0, y[1:m] >= 0, Int)

  ## Constraints

  ## Objective
  @objective(m0, Min, sum(y[i] for i in 1:m))

  #resolution
  optimize!(m0)

  y_star = JuMP.getvalue.( m0[:y] )

  valid = false

  while valid == false
    #### sous-probleme ####
    # Create the model
    m1 = Model(CPLEX.Optimizer)

    ## Variables
    @variable(m1, v[1:m+n] >= 0)

    ## Constraints
    @constraint(m1, [i in 1:m], - v[i] - v[adj[2*i-1]] + v[adj[2*i]] <= 0)
    @constraint(m1, [i in 1:m], - v[i] + v[adj[2*i-1]] - v[adj[2*i]] <= 0)
    @constraint(m1, v[1] == 0) # on saute la source
    @constraint(m1, sum(v[i] for i in 1:m+n) == 1)

    ## Objective
    @objective(m1, Max, - sum(b*y_star[ij]*v[n+ij] for ij in 1:m) + sum(demande[i]*v[i] for i in 1:n))

    #resolution
    optimize!(m1)

    v_star = JuMP.getvalue.( m1[:v] )

    if - sum(b*y_star[ij]*v_star[n+ij] for ij in 1:m) + sum(demande[i]*v_star[i] for i in 1:n) > 0 + 1e-8

      #### adding Constraints
      @constraint(m0, - sum(b*y[ij]*v_star[n+ij] for ij in 1:m) + sum(demande[i]*v_star[i] for i in 1:n) <= 0) # (4b)

      optimize!(m0)
      y_star = JuMP.getvalue.( m0[:y] )

    else
      valid = true
    end
  end
  println(y_star)
  println(JuMP.getobjectivevalue.( m0 ))
end

max_stars_repo_name: boumpteryx/MLA

max_issues_repo_name: boumpteryx/MLA

max_forks_repo_name: boumpteryx/MLA


text:
# This file was generated, do not modify it. # hide
nMissings = length(findall(x -> ismissing(x), data.commissioning_year))


text:
printstyled("  FastIO\n", color=:light_green)

# create testfile
testfile = "test.dat"
io = open(testfile, "w")
x = rand(UInt8, 1000)
x[255] = 0x0a
write(io, x)
close(io)

function fastio_test!(io::IO)
  T = [Bool, Int8, UInt8, Int16, UInt16, Float16, Int32, UInt32, Float32, Int64, UInt64, Float64]
  X = Array{Any, 1}(undef, length(T))
  Y = similar(X)
  r = rand(1:100)

  # =====================================================================
  # fastpos, fasteof
  @test fastpos(io) == position(io) == 0
  @test fasteof(io) == eof(io) == false
  seekend(io); p1=position(io); seekstart(io)
  fastseekend(io); p2=fastpos(io); seekstart(io)
  @test p1==p2

  # =====================================================================
  # test that a sequence of read, skip, and seek operations is identical
  skip(io, r)
  for (i,t) in enumerate(T)
    X[i] = read(io, t)
  end
  p1 = position(io)
  seek(io, r)
  p2 = position(io)
  seekstart(io)

  fastskip(io, r)
  for (i,t) in enumerate(T)
    Y[i] = fastread(io, t)
  end
  p3 = fastpos(io)
  fastseek(io, r)
  p4 = fastpos(io)

  @test p1==p3
  @test p2==p4
  for i=1:length(X)
    if isnan(X[i])
      @test isnan(Y[i])
    else
      @test X[i]==Y[i]
    end
  end

  seekstart(io)
  s1 = readline(io)
  seekstart(io)
  s2 = readline(io)
  @test s1 == s2
  return nothing
end

# =====================================================================
# done with file; read into buffer and repeat tests
printstyled("    IOStream\n", color=:light_green)
io = open(testfile, "r")
fastio_test!(io)
close(io)

printstyled("    generic IO\n", color=:light_green)
readbuf = read(testfile)
io = IOBuffer(readbuf)
fastio_test!(io)
close(io)

safe_rm(testfile)


text:
module EpidemicModel

    using Agents
    using DataFrames
    using LightGraphs
    using Random
    using StatsBase
    using Feather
    using JSON
    using Distributions
    using GraphIO
    using Distributed
    using Pipe
    using Query
    using Statistics
    using ProgressMeter

    include("citizen.jl")
    include("config.jl")
    include("initialization.jl")
    include("seir.jl")
    include("agent_step.jl")
    include("model_step.jl")
    include("formatting.jl")
    include("simulation.jl")
    include("io.jl")

    export Config
    export initialize_model
    export run_model!
    export run_experiment
    export write_dataframe
    export write_configs
    export write_networks
    export compress_data

end


text:
export getBoundModel,getBoundModelInv
function getBoundModel(m,a,b)
	d = (b.-a)./2.0;
	mid = (b.+a)./2.0;
	v = (1.0./d).*(m .- mid);
	s = d.*(tanh.(v).+1) .+ a;
	dsdm = sparse(Diagonal(sech.(v).^2));
	return s, dsdm
end

function getBoundModelInv(s,a,b)
	d = (b.-a)./2.0;
	mid = (b.+a)./2.0;
    m = atanh.((s.-a).*(1.0./d).-1) .* d .+ mid;
	return m
end


text:
# Archimedean copulas

"""
  arch_gen(copula::String, r::Matrix{Real}, θ::Real)

Auxiliary function used to generate data from archimedean copula (clayton, gumbel, frank or amh)
parametrized by a single parameter θ given a matrix of independent [0,1] distributerd
random vectors.

```jldoctest
  julia> arch_gen("clayton", [0.2 0.6 0.9; 0.4 0.5 0.8], 2.)
  2×2 Array{Real,2}:
   0.675778  0.851993
   0.687482  0.736394
```
"""

function arch_gen(copula, r, θ; rng)
  T = typeof(θ)
  t = size(r,1)
  if copula == "clayton"
    U = zeros(T, t, size(r,2)-1)
    for j in 1:t
      U[j,:] = clayton_gen(r[j,:], θ)
    end
    return U
  elseif copula == "amh"
    U = zeros(T ,t, size(r,2)-1)
    for j in 1:t
      U[j,:] = amh_gen(r[j,:], θ)
    end
    return U
  elseif copula == "frank"
    U = zeros(T, t, size(r,2)-1)
    w = logseriescdf(1-exp(-θ))
    for j in 1:t
      U[j,:] = frank_gen(r[j,:], θ, w)
    end
    return U
  else
    U = zeros(T, t, size(r,2)-1)
    u = r[:,end]
    p = invperm(sortperm(u))
    v = [levyel(θ, rand(rng), rand(rng)) for i in 1:length(u)]
    v = sort(v)[p]
    for j in 1:t
      U[j,:] = -log.(r[j,1:end-1])./v[j]
    end
    return exp.(-U.^(1/θ))
  end
end

"""
      gumbel_gen(r::Vector{Real}, θ::Real)

Given a vector of random numbers r of size n+2, return the sample from the Gumbel
copula parametrised by θ of the size n.
"""
function gumbel_gen(r, θ)
    u = -log.(r[1:end-2])./levyel(θ, r[end-1], r[end])
    return exp.(-u.^(1/θ))
end

"""
    clayton_gen(r::Vector{Real}, θ::Real)

Given a vector of random numbers r of size n+1, return the sample from the Clayton
copula parametrised by θ of the size n.
"""
function clayton_gen(r, θ)
    gamma_vec = gamma_inc_inv(1/θ, r[end], 1-r[end])
    u = -log.(r[1:end-1])./gamma_vec
    return (1 .+ u).^(-1/θ)
end

"""
    amh_gen(r::Vector{Real}, θ::Real)

Given a vector of random numbers r of size n+1, return the sample from the AMH
copula parametrised by θ of the size n.
"""
function amh_gen(r, θ)
    u = -log.(r[1:end-1])./(1 .+quantile.(Geometric(1-θ), r[end]))
    return (1-θ) ./(exp.(u) .-θ)
end

"""
    frank_gen(r::Vector{Real}, θ::Real, logseries::Vector{Real})

Given a vector of random numbers r of size n+1, return the sample from the Frank
copula parametrised by θ of the size n. Axiliary logseries is a vector of the
logseries sequence priorly computed.
"""
function frank_gen(r, θ, logseries)
    logseriesquantile = findlast(logseries .< r[end])
    u = -log.(r[1:end-1])./logseriesquantile
    return -log.(1 .+exp.(-u) .*(exp(-θ)-1)) ./θ
end

"""
    Gumbel_cop

Fields:
  - n::Int - number of marginals
  - θ::Real - parameter

Constructor

        Gumbel_cop(n::Int, θ::Real)

The Gumbel n variate copula is parameterized by θ::Real ∈ [1, ∞), supported for n::Int ≧ 2.

Constructor

    Gumbel_cop(n::Int, θ::Real, cor::Type{<:CorrelationType})

For computing copula parameter from expected correlation use empty type cor::Type{<:CorrelationType} where
SpearmanCorrelation <:CorrelationType and KendallCorrelation<:CorrelationType. If used cor put expected correlation in the place of θ  in the constructor.
The copula parameter will be computed then. The correlation must be greater than zero.

```jldoctest

julia> Gumbel_cop(4, 3.)
Gumbel_cop(4, 3.0)

julia> Gumbel_cop(4, .75, KendallCorrelation)
Gumbel_cop(4, 4.0)
```
"""
struct Gumbel_cop{T} <: Copula{T}
  n::Int
  θ::T
  function(::Type{Gumbel_cop})(n::Int, θ::T) where T <: Real
      n >= 2 || throw(DomainError("not supported for n < 2"))
      testθ(θ, "gumbel")
      new{T}(n, θ)
  end
  function(::Type{Gumbel_cop})(n::Int, ρ::T, cor::Type{<:CorrelationType}) where T <: Real
      n >= 2 || throw(DomainError("not supported for n < 2"))
      θ = getθ4arch(ρ, "gumbel", cor)
      new{T}(n, θ)
  end
end

"""
    simulate_copula!(U::Matrix{Real}, copula::Gumbel_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Given the preallocated output U, Returns size(U,1) realizations from the Gumbel copula -  Gumbel_cop(n, θ)
N.o. marginals is size(U,2), requires size(U,2) == copula.n

```jldoctest
julia> Random.seed!(43);

julia> U = zeros(2,3)
2×3 Array{Float64,2}:
 0.0  0.0  0.0
 0.0  0.0  0.0

julia> simulate_copula!(U, Gumbel_cop(3, 1.5))

julia> U
2×3 Array{Float64,2}:
 0.740038  0.918928  0.950674
 0.637826  0.483514  0.123949
```
"""
function simulate_copula!(U, copula::Gumbel_cop{T}; rng = Random.GLOBAL_RNG) where T
    θ = copula.θ
    n = copula.n
    size(U, 2) == n || throw(AssertionError("n.o. margins in pre allocated output and copula not equal"))
    for j in 1:size(U,1)
      u = rand(rng, T, n+2)
      U[j,:] = gumbel_gen(u, θ)
    end
end

"""
    Gumbel_cop_rev

Fields:
  - n::Int - number of marginals
  - θ::Real - parameter

Constructor

    Gumbel_cop_rev(n::Int, θ::Real)

The reversed Gumbel copula (reversed means u → 1 .- u),
parameterized by θ::Real ∈ [1, ∞), supported for n::Int ≧ 2.

Constructor

    Gumbel_cop_rev(n::Int, θ::Real, cor::Type{<:CorrelationType})

For computing copula parameter from expected correlation use empty type cor::Type{<:CorrelationType} where
SpearmanCorrelation <:CorrelationType and KendallCorrelation<:CorrelationType. If used cor put expected correlation in the place of θ  in the constructor.
The copula parameter will be computed then. The correlation must be greater than zero.

```jldoctest
julia> c = Gumbel_cop_rev(4, .75, KendallCorrelation)
Gumbel_cop_rev(4, 4.0)

julia> Random.seed!(43);

julia> simulate_copula(2, c)
2×4 Array{Float64,2}:
 0.963524   0.872108  0.816626  0.783637
 0.0954475  0.138451  0.13593   0.0678172
```
"""
struct Gumbel_cop_rev{T} <: Copula{T}
  n::Int
  θ::T
  function(::Type{Gumbel_cop_rev})(n::Int, θ::T) where T <: Real
      n >= 2 || throw(DomainError("not supported for n < 2"))
      testθ(θ, "gumbel")
      new{T}(n, θ)
  end
  function(::Type{Gumbel_cop_rev})(n::Int, ρ::T, cor::Type{<:CorrelationType}) where T <: Real
      n >= 2 || throw(DomainError("not supported for n < 2"))
      θ = getθ4arch(ρ, "gumbel", cor)
      new{T}(n, θ)
  end
end

"""
    simulate_copula!(U::Matrix{Real}, copula::Gumbel_cop_rev; rng::AbstractRNG = Random.GLOBAL_RNG)

Given the preallocated output U, Returns size(U,1) realizations from the reversed Gumbel copula -  Gumbel_cop_rev(n, θ)
N.o. marginals is size(U,2), requires size(U,2) == copula.n

```jldoctest
julia> Random.seed!(43);

julia> U = zeros(2,3)
2×3 Array{Float64,2}:
 0.0  0.0  0.0
 0.0  0.0  0.0

julia> simulate_copula!(U, Gumbel_cop_rev(3, 1.5))

julia> U
2×3 Array{Flaot64,2}:
 0.259962  0.081072  0.0493259
 0.362174  0.516486  0.876051
```
"""
function simulate_copula!(U, copula::Gumbel_cop_rev{T}; rng = Random.GLOBAL_RNG) where T
    θ = copula.θ
    n = copula.n
    size(U, 2) == n || throw(AssertionError("n.o. margins in pre allocated output and copula not equal"))
    for j in 1:size(U,1)
      u = rand(rng, T, n+2)
      U[j,:] = 1 .- gumbel_gen(u, θ)
    end
end

"""
    Clayton_cop

Fields:
  - n::Int - number of marginals
  - θ::Real - parameter

Constructor

    Clayton_cop(n::Int, θ::Real)

The Clayton n variate copula parameterized by θ::Real, such that θ ∈ (0, ∞) for n > 2 and θ ∈ [-1, 0) ∪ (0, ∞) for n = 2,
supported for n::Int ≥ 2.

Constructor

    Clayton_cop(n::Int, θ::Real, cor::Type{<:CorrelationType})

For computing copula parameter from expected correlation use empty type cor::Type{<:CorrelationType} where
SpearmanCorrelation <:CorrelationType and KendallCorrelation<:CorrelationType. If used cor put expected correlation in the place of θ  in the constructor.
The copula parameter will be computed then. The correlation must be greater than zero.

```jldoctest
julia> Clayton_cop(4, 3.)
Clayton_cop(4, 3.0)

julia> Clayton_cop(4, 0.9, SpearmanCorrelation)
Clayton_cop(4, 5.5595567742323775)
```
"""
struct Clayton_cop{T} <: Copula{T}
  n::Int
  θ::T
  function(::Type{Clayton_cop})(n::Int, θ::T) where T <: Real
      n >= 2 || throw(DomainError("not supported for n < 2"))
      if n > 2
        testθ(θ, "clayton")
      else
        (θ >= -1.) & (θ != 0.) || throw(DomainError("bivariate Clayton not supported for θ < -1 or θ = 0"))
      end
      new{T}(n, θ)
  end
  function(::Type{Clayton_cop})(n::Int, ρ::T, cor::Type{<:CorrelationType}) where T <: Real
      n >= 2 || throw(DomainError("not supported for n < 2"))
      θ = getθ4arch(ρ, "clayton", cor)
      new{T}(n, θ)
  end
end

"""
    simulate_copula!(U::Matrix{Real}, copula::Clayton_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Given the preallocated output U, Returns t realizations from the Clayton copula - Clayton_cop(n, θ)
N.o. marginals is size(U,2), requires size(U,2) == copula.n.
N.o. realisations is size(U,1).

```jldoctest
julia> Random.seed!(43);

julia> U = zeros(3,2)
3×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0
 0.0  0.0

julia> simulate_copula!(U, Clayton_cop(2, 1.))

julia> U
3×2 Array{Float64,2}:
 0.562482  0.896247
 0.968953  0.731239
 0.749178  0.38015

julia> U = zeros(2,2)
2×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0

julia> Random.seed!(43);

julia> simulate_copula!(U, Clayton_cop(2, -.5))

julia> U
2×2 Array{Float64,2}:
 0.180975  0.818017
 0.888934  0.863358
```
"""
function simulate_copula!(U, copula::Clayton_cop{T}; rng = Random.GLOBAL_RNG) where T
    θ = copula.θ
    n = copula.n
    size(U, 2) == n || throw(AssertionError("n.o. margins in pre allocated output and copula not equal"))
    if (n == 2) & (θ < 0)
        simulate_copula!(U, Chain_of_Archimedeans([θ], "clayton"); rng = rng)
    else
        for j in 1:size(U,1)
            u = rand(rng, T, n+1)
            U[j,:] = clayton_gen(u, θ)
        end
    end
end

"""
    Clayton_cop_rev

Fields:
- n::Int - number of marginals
- θ::Real - parameter

Constructor

    Clayton_cop_rev(n::Int, θ::Real)

The reversed Clayton copula parameterized by θ::Real (reversed means u → 1 .- u).
Domain: θ ∈ (0, ∞) for n > 2 and θ ∈ [-1, 0) ∪ (0, ∞) for n = 2,
supported for n::Int ≧ 2.

Constructor

    Clayton_cop_rev(n::Int, θ::Real, cor::Type{<:CorrelationType})

For computing copula parameter from expected correlation use empty type cor::Type{<:CorrelationType} where
SpearmanCorrelation <:CorrelationType and KendallCorrelation<:CorrelationType. If used cor put expected correlation in the place of θ  in the constructor.
The copula parameter will be computed then. The correlation must be greater than zero.

```jldoctest

julia> Clayton_cop_rev(4, 3.)
Clayton_cop_rev(4, 3.0)

julia> Clayton_cop_rev(4, 0.9, SpearmanCorrelation)
Clayton_cop_rev(4, 5.5595567742323775)

```
"""
struct Clayton_cop_rev{T} <: Copula{T}
  n::Int
  θ::T
  function(::Type{Clayton_cop_rev})(n::Int, θ::T) where T <: Real
      n >= 2 || throw(DomainError("not supported for n < 2"))
      if n > 2
        testθ(θ, "clayton")
      else
        (θ >= -1.) & (θ != 0.) || throw(DomainError("bivariate Clayton not supported for θ < -1 or θ = 0"))
      end
      new{T}(n, θ)
  end
  function(::Type{Clayton_cop_rev})(n::Int, ρ::T, cor::Type{<:CorrelationType}) where T <: Real
      n >= 2 || throw(DomainError("not supported for n < 2"))
      θ = getθ4arch(ρ, "clayton", cor)
      new{T}(n, θ)
  end
end

"""
    simulate_copula!(U::Matrix{Real}, copula::Clayton_cop_rev; rng::AbstractRNG = Random.GLOBAL_RNG)

Given the preallocated output U, Returns size(U,1) realizations from the reversed Clayton copula - Clayton_cop_rev(n, θ)
N.o. marginals is size(U,2), requires size(U,2) == copula.n

```jldoctest
julia> Random.seed!(43);

julia> U = zeros(2,2)
2×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0

julia> simulate_copula!(U, Clayton_cop_rev(2, -0.5))

julia> U
2×2 Array{Float64,2}:
 0.819025  0.181983
 0.111066  0.136642

 julia> Random.seed!(43);

julia> U = zeros(2,2)
2×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0

julia> simulate_copula!(U, Clayton_cop_rev(2, 2.))

julia> U
2×2 Array{Float64,2}:
 0.347188   0.087281
 0.0257036  0.212676
```
"""
function simulate_copula!(U, copula::Clayton_cop_rev{T}; rng = Random.GLOBAL_RNG) where T
    θ = copula.θ
    n = copula.n
    size(U, 2) == n || throw(AssertionError("n.o. margins in pre allocated output and copula not equal"))
    if (n == 2) & (θ < 0)
        for j in 1:size(U,1)
          u = rand(rng, T)
          w = rand(rng, T)
          U[j,:] = 1 .- hcat(u, rand2cop(u, θ, "clayton", w))
        end
    else
        for j in 1:size(U,1)
            u = rand(rng, T, n+1)
            U[j,:] = 1 .- clayton_gen(u, θ)
        end
    end
end

"""
    AMH_cop

Fields:
- n::Int - number of marginals
- θ::Real - parameter

Constructor

    AMH_cop(n::Int, θ::Real)

The Ali-Mikhail-Haq copula parameterized by θ, domain: θ ∈ (0, 1) for n > 2 and  θ ∈ [-1, 1] for n = 2.

Constructor

    AMH_cop(n::Int, θ::Real, cor::Type{<:CorrelationType})

For computing copula parameter from expected correlation use empty type cor::Type{<:CorrelationType} where
SpearmanCorrelation <:CorrelationType and KendallCorrelation<:CorrelationType. If used cor put expected correlation in the place of θ  in the constructor.
The copula parameter will be computed then. The correlation must be greater than zero.
Such correlation must be grater than zero and limited from above due to the θ domain.
            - Spearman correlation must be in range (0, 0.5)
            - Kendall correlation must be in range (0, 1/3)

```jldoctest
julia> AMH_cop(4, .3)
AMH_cop(4, 0.3)

julia> AMH_cop(4, .3, KendallCorrelation)
AMH_cop(4, 0.9999)

```
"""
struct AMH_cop{T} <: Copula{T}
  n::Int
  θ::T
  function(::Type{AMH_cop})(n::Int, θ::T) where T <: Real
      n >= 2 || throw(DomainError("not supported for n < 2"))
      if n > 2
        testθ(θ, "amh")
      else
      1 >=  θ >= -1 || throw(DomainError("bivariate AMH not supported for θ > 1 or θ < -1"))
      end
      new{T}(n, θ)
  end
  function(::Type{AMH_cop})(n::Int, ρ::T, cor::Type{<:CorrelationType}) where T <: Real
      n >= 2 || throw(DomainError("not supported for n < 2"))
      θ = getθ4arch(ρ, "amh", cor)
      new{T}(n, θ)
  end
end

"""
    simulate_copula!(U::Matrix{Real}, copula::AMH_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Given the preallocated output U, Returns size(U,1) realizations from the Ali-Mikhail-Haq copula- AMH_cop(n, θ)
N.o. marginals is size(U,2), requires size(U,2) == copula.n

```jldoctest
julia> Random.seed!(43);

julia> simulate_copula!(U, AMH_cop(2, -0.5))

julia> U
4×2 Array{Float64,2}:
 0.180975  0.820073
 0.888934  0.886169
 0.408278  0.919572
 0.828727  0.335864
```
"""
function simulate_copula!(U, copula::AMH_cop{T}; rng = Random.GLOBAL_RNG) where T
  n = copula.n
  θ = copula.θ
  size(U, 2) == n || throw(AssertionError("n.o. margins in pre allocated output and copula not equal"))
  if (θ in [0,1]) | (n == 2)*(θ < 0)
      simulate_copula!(U, Chain_of_Archimedeans([θ], "amh"); rng = rng)
  else
    for j in 1:size(U,1)
      u = rand(rng, T, n+1)
      U[j,:] = amh_gen(u, θ)
    end
  end
end

"""
    AMH_cop_rev

Fields:
  - n::Int - number of marginals
  - θ::Real - parameter

Constructor

    AMH_cop_rev(n::Int, θ::Real)

The reversed Ali-Mikhail-Haq copula parametrized by θ, i.e.
such that the output is 1 .- u, where u is modelled by the corresponding AMH copula.
Domain: θ ∈ (0, 1) for n > 2 and  θ ∈ [-1, 1] for n = 2.

Constructor

    AMH_cop_rev(n::Int, θ::Real, cor::Type{<:CorrelationType})

For computing copula parameter from expected correlation use empty type cor::Type{<:CorrelationType} where
SpearmanCorrelation <:CorrelationType and KendallCorrelation<:CorrelationType. If used cor put expected correlation in the place of θ  in the constructor.
The copula parameter will be computed then. The correlation must be greater than zero.
Such correlation must be grater than zero and limited from above due to the θ domain.
              - Spearman correlation must be in range (0, 0.5)
              -  Kendall correlation must be in range (0, 1/3)

```jldoctest
julia> AMH_cop_rev(4, .3)
AMH_cop_rev(4, 0.3)
```
"""
struct AMH_cop_rev{T} <: Copula{T}
  n::Int
  θ::T
  function(::Type{AMH_cop_rev})(n::Int, θ::T) where T <: Real
      n >= 2 || throw(DomainError("not supported for n < 2"))
      if n > 2
        testθ(θ, "amh")
      else
        1 >=  θ >= -1 || throw(DomainError("bivariate AMH not supported for θ > 1 or θ < -1"))
      end
      new{T}(n, θ)
  end
  function(::Type{AMH_cop_rev})(n::Int, ρ::T, cor::Type{<:CorrelationType}) where T <: Real
      n >= 2 || throw(DomainError("not supported for n < 2"))
      θ = getθ4arch(ρ, "amh", cor)
      new{T}(n, θ)
  end
end

"""
    simulate_copula!(U::Matrix{Real}, copula::AMH_cop_rev; rng::AbstractRNG = Random.GLOBAL_RNG)

Given the preallocated output U, Returns size(U,1) realizations from the reversed Ali-Mikhail-Haq copula a - AMH_cop_rev(n, θ)
N.o. marginals is size(U,2), requires size(U,2) == copula.n

```jldoctest
julia> Random.seed!(43);

julia> U = zeros(4,2)
4×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0
 0.0  0.0
 0.0  0.0

julia> simulate_copula!(U, AMH_cop_rev(2, 0.5))

julia> U
4×2 Array{Float64,2}:
 0.516061   0.116089
 0.0379356  0.334231
 0.292457   0.74958
 0.0845089  0.505477
```
"""
function simulate_copula!(U, copula::AMH_cop_rev{T}; rng = Random.GLOBAL_RNG) where T
    θ = copula.θ
    n = copula.n
    size(U, 2) == n || throw(AssertionError("n.o. margins in pre allocated output and copula not equal"))
    if (n == 2) & (θ < 0)
        for j in 1:size(U,1)
          u = rand(rng, T)
          w = rand(rng, T)
          U[j,:] = 1 .- hcat(u, rand2cop(u, θ, "amh", w))
        end
    else
        for j in 1:size(U,1)
            u = rand(rng, T, n+1)
            U[j,:] = 1 .- amh_gen(u, θ)
        end
    end
end

"""
    Frank_cop

Fields:
- n::Int - number of marginals
- θ::Real - parameter

Constructor

    Frank_cop(n::Int, θ::Real)

The Frank n variate copula parameterized by θ::Real.
Domain: θ ∈ (0, ∞) for n > 2 and θ ∈ (-∞, 0) ∪ (0, ∞) for n = 2,
supported for n::Int ≧ 2.

Constructor

    Frank_cop(n::Int, θ::Real, cor::Type{<:CorrelationType})

For computing copula parameter from expected correlation use empty type cor::Type{<:CorrelationType} where
SpearmanCorrelation <:CorrelationType and KendallCorrelation<:CorrelationType. If used cor put expected correlation in the place of θ  in the constructor.
The copula parameter will be computed then. The correlation must be greater than zero.

```jldoctest
julia> Frank_cop(2, -5.)
Frank_cop(2, -5.0)

julia> Frank_cop(4, .3)
Frank_cop(4, 0.3)
```
"""
struct Frank_cop{T} <: Copula{T}
  n::Int
  θ::T
  function(::Type{Frank_cop})(n::Int, θ::T) where T <: Real
      n >= 2 || throw(DomainError("not supported for n < 2"))
      if n > 2
        testθ(θ, "frank")
      else
        θ != 0 || throw(DomainError("bivariate frank not supported for θ = 0"))
      end
      new{T}(n, θ)
  end
  function(::Type{Frank_cop})(n::Int, ρ::T, cor::Type{<:CorrelationType}) where T <: Real
      n >= 2 || throw(DomainError("not supported for n < 2"))
      θ = getθ4arch(ρ, "frank", cor)
      new{T}(n, θ)
  end
end

"""
    simulate_copula!(U::Matrix{Real}, copula::Frank_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Given the preallocated output U, Returns size(U,1) realizations from the Frank copula- Frank_cop(n, θ)
N.o. marginals is size(U,2), requires size(U,2) == copula.n

```jldoctest
julia> U = zeros(4,2)
4×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0
 0.0  0.0
 0.0  0.0

julia> Random.seed!(43);

julia> simulate_copula!(U, Frank_cop(2, 3.5))

julia> U
4×2 Array{Float64,2}:
 0.650276  0.910212
 0.973726  0.789701
 0.690966  0.358523
 0.747862  0.29333
```
"""
function simulate_copula!(U, copula::Frank_cop{T}; rng = Random.GLOBAL_RNG) where T
  n = copula.n
  θ = copula.θ
  size(U, 2) == n || throw(AssertionError("n.o. margins in pre allocated output and copula not equal"))
  if (n == 2) & (θ < 0)
      simulate_copula!(U, Chain_of_Archimedeans([θ], "frank"); rng = rng)
  else
    w = logseriescdf(1-exp(-θ))
    for j in 1:size(U,1)
      u = rand(rng, T, n+1)
      U[j,:] = frank_gen(u, θ, w)
    end
  end
end

"""
  function testθ(θ::Real, copula::String)

Tests the parameter θ value for archimedean copula, returns void
"""
function testθ(θ, copula) where T
  if copula == "gumbel"
    θ >= 1 || throw(DomainError("gumbel copula not supported for θ < 1"))
  elseif copula == "amh"
    1 > θ > 0 || throw(DomainError("amh multiv. copula supported only for 0 < θ < 1"))
  else
    θ > 0 || throw(DomainError("generaton not supported for θ ≤ 0"))
  end
  Nothing
end

"""
  useρ(ρ::Real, copula::String)

Tests the available Spearman correlation for the Archimedean copula.

Returns Real, the copula parameter θ with the Spearman correlation ρ.
```jldoctest

julia> useρ(0.75, "gumbel")
2.294053859606698
```

"""
function useρ(ρ, copula)
  0 < ρ < 1 || throw(DomainError("Spearman correlation coeficiant must fulfill 0 < ρ < 1"))
  if copula == "amh"
    0 < ρ < 0.5 || throw(DomainError("Spearman correlation coeficiant must fulfill 0 < ρ < 0.5"))
  end
  ρ2θ(ρ, copula)
end

"""
  useτ(ρ::Real, copula::String)

Tests the available kendall's correlation for archimedean copula, returns Float,
corresponding copula parameter θ.

```jldoctest

julia> useτ(0.5, "clayton")
2.0
```
"""
function useτ(τ, copula)
  0 < τ < 1 || throw(DomainError("Kendall correlation coeficiant must fulfill 0 < τ < 1"))
  if copula == "amh"
    0 < τ < 1/3 || throw(DomainError("Kendall correlation coeficiant must fulfill 0 < τ < 1/3"))
  end
  τ2θ(τ, copula)
end

"""

    getθ4archgetθ4arch(ρ::Real, copula::String, cor::Type{SpearmanCorrelation})

    getθ4archgetθ4arch(ρ::Real, copula::String, cor::Type{KendallCorrelation})

Compute the copula parameter given the correlation, test if the parameter is in range.
Following types are supported: SpearmanCorrelation, KendallCorrelation

```jldoctest

julia> getθ4arch(0.5, "gumbel", SpearmanCorrelation)
1.541070420842913

julia> getθ4arch(0.5, "gumbel", KendallCorrelation)
2.0
```
"""
getθ4arch(ρ, copula, cor::Type{SpearmanCorrelation}) = useρ(ρ , copula)

getθ4arch(ρ, copula, cor::Type{KendallCorrelation}) = useτ(ρ , copula)


text:
@with_kw mutable struct WaterSource <: AgComponent
    name::String
    cost_per_ML::Union{Int64, Float64, AgParameter}
    cost_per_ha::Union{Int64, Float64, AgParameter}
    yearly_cost::Union{Int64, Float64, AgParameter}
    pump::Pump
    head::Union{Int64, Float64, AgParameter}
    allocation::Union{Int64, Float64, AgParameter}
    entitlement::Union{Int64, Float64, AgParameter}
end

function pump_cost_per_ML(ws::WaterSource, flow_rate_Lps::Float64)::Float64
    return pumping_costs_per_ML(ws.pump, flow_rate_Lps, ws.head)
end

function usage_costs(ws::WaterSource, water_used_ML::Float64)::Float64
    return ws.cost_per_ML * water_used_ML
end

function subtotal_costs(ws::WaterSource, area::Float64, water_used_ML::Float64)::Float64
    usage_fee::Float64 = usage_costs(ws, water_used_ML)
    area_fee::Float64 = ws.cost_per_ha * area
    return ws.yearly_cost + usage_fee + area_fee
end

function create(ws_specs::Dict{Symbol, Dict},
                pump_specs::Dict{Symbol, Dict};
                id_prefix::Union{String, Nothing}=nothing)::Array{WaterSource}
    ws::Array{WaterSource} = []
    for (k, v) in ws_specs
        pump_name = v[:name]  # pump will have same name as the water_source
        v[:pump] = create(pump_specs[Symbol(pump_name)])

        ws_i::WaterSource = create(v)
        push!(ws, ws_i)
    end

    return ws
end

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
import MPI

precompile(Tuple{typeof(Diffusion.julia_diffusion), Ptr{Floa64}, Cint, Cint, Cint, MPI.MPI_Comm})

max_stars_repo_name: JuliaLabs/CallingJuliaFromCWithMPI

max_issues_repo_name: JuliaLabs/CallingJuliaFromCWithMPI

max_forks_repo_name: JuliaLabs/CallingJuliaFromCWithMPI


text:
export isless

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

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

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

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

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

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

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

    return( true )
end


text:
module DemoPackageCT

# Write your package code here.
print_greeting() = printstyled("Hello, world!"; color=:magenta)

end


text:
mutable struct CustomModel <: AbstractModel
    settings::Settings
    constraints::Vector{Constraint}
    obj::CustomObjective
    output::Output
    CustomModel(settings, constraints, obj, output) =
        new(settings, constraints, obj, output)
    CustomModel() = new(Settings(), Constraint[], CustomObjective(), Output())
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule jemalloc_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("jemalloc")
JLLWrappers.@generate_main_file("jemalloc", UUID("454a8cc1-5e0e-5123-92d5-09b094f0e876"))
end  # module jemalloc_jll


text:
module StructViews

include("abstractview.jl")
include("structview.jl")
include("fieldview.jl")

export StructView
export FieldView

end # module


text:
#=
  From
  "Probabilistic logic programming and its applications"
  Luc De Raedt, Leuven
  https://www.youtube.com/watch?v=3lnVBqxjC88
  @ 3:44

  """
  Mike has a bag of marbles with 4 white, 8 blue, and
  6 red marbles. He pulls out one marble from the bag
  and it is red. What is the probability that the
  second marble he pulls out of the bag is white?

  The answer is 0.234931.
  """

   draw1: The probability of drawing white is the probability for 1.0
   Distributions of variable draw[1]
   red        =>  100000  (1.000000)

   Distributions of variable draw[2]
   blue       =>   47339  (0.473390)
   red        =>   29353  (0.293530)
   white      =>   23308  (0.233080)

  Cf ~/cplint/bag_of_marbles.pl
     ~/blog/bag_of_marbles.blog
     ~/psi/bag_of_marbles.psi
     ~/webppl/bag_of_marbles.wppl

  Cf bag_of_marbles.jl. This current model use a for loop.

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

#=
draw1: The probability of drawing white is the probability for 1.0

Distributions of variable draw[2] (num:0)
2.00000 =>  188281  (0.470703)
3.00000 =>  117185  (0.292963)
1.00000 =>   94534  (0.236335)
=#

@model function bag_of_marbles()

    white = 1
    blue  = 2
    red   = 3

    start = [4,8,6]

    draw = tzeros(2)
    for t in 1:2
        # We have to use sum(Int8[...]) since sum([]) throws an error
        draw[t] ~ Categorical(
                     simplex([
                            start[white] - sum([draw[i] == white ? 1 : 0 for i in 1:t-1]),
                            start[blue]  - sum([draw[i] == blue  ? 1 : 0 for i in 1:t-1]),
                            start[red]   - sum([draw[i] == red   ? 1 : 0 for i in 1:t-1])
                            ])
                            )
    end

    true ~ Dirac(draw[1] == red)

end

model = bag_of_marbles()
num_chns = 4

# chns = sample(model, Prior(), 100_000)
chns = sample(model, MH(), 100_000)
# chns = sample(model, PG(15), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(),10_000)

display(chns)

println("\ndraw1: The probability of drawing white is the probability for 1.0")
show_var_dist_pct(chns,Symbol("draw[1]"),["white","blue","red"])
show_var_dist_pct(chns,Symbol("draw[2]"),["white","blue","red"])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
test_render_figures(format, chunk) = Weave.render_figures(get_format(format), chunk)

# Make a dummy codehunk with figure
chunk = Weave.CodeChunk("plot(x)", 1, 1, "", Dict())
options = merge(Weave.get_chunk_defaults(), chunk.options)
merge!(chunk.options, options)
chunk.figures = ["figs/figures_plot1.png"]

@test test_render_figures("md2tex", chunk) == "\\includegraphics{figs/figures_plot1.png}\n"
@test test_render_figures("texminted", chunk) == "\\includegraphics{figs/figures_plot1.png}\n"
@test test_render_figures("pandoc", chunk) == "![](figs/figures_plot1.png)\\ \n\n"
@test test_render_figures("github", chunk) == "![](figs/figures_plot1.png)\n"
@test test_render_figures("hugo", chunk) == "{{
}}"
@test test_render_figures("multimarkdown", chunk) == "![][figs/figures_plot1.png]\n\n[figs/figures_plot1.png]: figs/figures_plot1.png \n"
@test test_render_figures("md2html", chunk) == "\n"

chunk.options[:out_width] = "100%"
@test test_render_figures("asciidoc", chunk) == "image::figs/figures_plot1.png[width=100%]\n"
@test test_render_figures("rst", chunk) == ".. image:: figs/figures_plot1.png\n   :width: 100%\n\n"

chunk.options[:fig_cap] = "Nice plot"
@test test_render_figures("pandoc", chunk) == "![Nice plot](figs/figures_plot1.png){width=100%}\n"
@test test_render_figures("md2tex", chunk) == "\\begin{figure}[!h]\n\\center\n\\includegraphics[width=1.0\\linewidth]{figs/figures_plot1.png}\n\\caption{Nice plot}\n\\end{figure}\n"
@test test_render_figures("md2html", chunk) == "
\n\n
Nice plot
\n
\n"
@test test_render_figures("rst", chunk) == ".. figure:: figs/figures_plot1.png\n   :width: 100%\n\n   Nice plot\n\n"
@test test_render_figures("multimarkdown", chunk) == "![Nice plot][figs/figures_plot1.png]\n\n[figs/figures_plot1.png]: figs/figures_plot1.png width=100%\n"
@test test_render_figures("asciidoc", chunk) == "image::figs/figures_plot1.png[width=100%,title=\"Nice plot\"]"

chunk.options[:label] = "somefig"
@test test_render_figures("pandoc", chunk) == "![Nice plot](figs/figures_plot1.png){width=100% #fig:somefig}\n"


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1TopologySelectorTerm <: SwaggerModel
    matchLabelExpressions::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1TopologySelectorLabelRequirement} } # spec name: matchLabelExpressions

    function IoK8sApiCoreV1TopologySelectorTerm(;matchLabelExpressions=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1TopologySelectorTerm, Symbol("matchLabelExpressions"), matchLabelExpressions)
        setfield!(o, Symbol("matchLabelExpressions"), matchLabelExpressions)
        o
    end
end # type IoK8sApiCoreV1TopologySelectorTerm

const _property_map_IoK8sApiCoreV1TopologySelectorTerm = Dict{Symbol,Symbol}(Symbol("matchLabelExpressions")=>Symbol("matchLabelExpressions"))
const _property_types_IoK8sApiCoreV1TopologySelectorTerm = Dict{Symbol,String}(Symbol("matchLabelExpressions")=>"Vector{IoK8sApiCoreV1TopologySelectorLabelRequirement}")
Base.propertynames(::Type{ IoK8sApiCoreV1TopologySelectorTerm }) = collect(keys(_property_map_IoK8sApiCoreV1TopologySelectorTerm))
Swagger.property_type(::Type{ IoK8sApiCoreV1TopologySelectorTerm }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1TopologySelectorTerm[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1TopologySelectorTerm }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1TopologySelectorTerm[property_name]

function check_required(o::IoK8sApiCoreV1TopologySelectorTerm)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1TopologySelectorTerm }, name::Symbol, val)
end


text:
"""
    zero_one_pack!(capacity::N, weights::V, values::V, dp::V) where {N <: Number,V <: AbstractVector}
This does 0-1 (each item can be chosen only once) knapsack :
pack capacity = capacity
weight of each item = weights
value of each item = values
dp array is what the function works on
It returns the ans (dp[capacity])
```
julia> dp=zeros(Int,30)
julia> zero_one_pack!(20,[1,3,11],[2,5,30],dp)
37
```
"""
function dp_zero_one_pack!(capacity::N, weights::V, values::V, dp::V
) where {N <: Number,V <: AbstractVector}
    for i in 1:length(weights)
        j = capacity
        while j > weights[i] # reversed loop
            dp[j] = max(dp[j], dp[j - weights[i]] + values[i])
            j -= 1
        end
        dp[weights[i]] = max(dp[weights[i]], values[i]) # dp[j]=max(dp[j],dp[0]+values[i])
    end
    return dp[capacity]
end

"""
This does complete/infinite (each item can be chosen infinite times) knapsack :
pack capacity = capacity
weight of each item = weights
value of each item = values
dp array is what the function works on
It returns the ans (dp[capacity])
```
julia> dp=zeros(Int,30)
julia> complete_pack!(20,[1,2,9],[1,3,20],dp)
43
```
"""
function dp_complete_pack!(capacity::N, weights::V, values::V, dp::V
) where {N <: Number,V <: AbstractVector}
    for i in 1:length(weights)
        dp[weights[i]] = max(dp[weights[i]], values[i]) # dp[j]=max(dp[j],dp[0]+values[i])
        for j in weights[i] + 1:capacity
            dp[j] = max(dp[j], dp[j - weights[i]] + values[i])
        end
    end
    return dp[capacity]
end


text:
# basisfunctions.jl

"""
    fourier(k::Int)

Implements the Fourier basis of functions ϕ_k, defined, if k is odd, by
ϕ_k(x)=sqrt(2)sin((k+1)π*x) and if k is even, by ϕ_k(x)=sqrt(2)cos(k*x).

# Examples

```julia
using PyPlot
clf()
x = -1.0:0.01:2.0
y = fourier(3).(x)
plot(x,y)
```
"""
function fourier(k::Int)
  k ≥ 1 || error("k should be positive")
  sqrttwo = sqrt(2.0)
  if k % 2 == 1
    return x -> sqrttwo * sinpi((k+1)*x)
  else
    return  x -> sqrttwo * cospi(k*x)
  end
end

"""
    faberschauderone(x::Float64)

Implements the first Faber-Schauder function defined by 1-2x for 0≤x≤1/2 and
-1+2x for 0.5≤x≤1, and is 1-periodic extended to all x∈R.

See also: faberschauder

#Example
```julia
x=-2.:0.001:2
y=faberschauderone.(x)
using PyPlot
clf()
plot(x,y)
```
"""
function faberschauderone(x::Float64)
    y = mod(x, 1.0)
    if 0≤y≤0.5
        return 1.0 - 2*y
    else
        return -1.0 + 2*y
    end
end

"""
    faberschauder(j::Int, k::Int)

Implements the k-th Faber-Schauder function of level j. Here, j≥0 and 1≤k≤2^j.
It is a one periodic function and defined on [0,1] by 2^(j+1)(x-(k-1)2^j) on
(k-1)2^j≤x≤(2k-1)2^(j+1) and 1 - 2^(j+1)(x-(2k-1)2^(j+1)) on
[(2k-1)2^(j+1), k2^j] and zero outside these intervals.

#Example
```julia
j=3
x=-2.:0.001:2.
using PyPlot
clf()
for k in 1:2^j
    y = faberschauder(j,k).(x)
    plot(x,y)
end
```
"""
function faberschauder(j::Int,k::Int)
  j ≥ 0 || error("j should be a nonnegative integer.")
  1≤k≤2^j || error("k should be an integer between 1 and 2^j.")
  j_float = Float64(j)
  return function(x::Float64)
    y = mod(x,1.0)
    if y ≤ (k-1)*2.0^(-j_float) || y ≥ k*2.0^(-j_float)
      return 0.0
    elseif ((k-1)*2.0^(-j_float) < y ≤ (2k-1)*2.0^(-j_float-1.0))
      return 2^(j+1)*(y-(k-1)*2.0^(-j_float))
    else
      return 1.0 - 2^(j+1)*(y-(2k-1)*2.0^(-j_float-1.0))
    end
  end
end


text:
using DiffEqSensitivity, Flux, Zygote, StochasticDiffEq, Test

function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end
function lotka_volterra(u,p,t)
  x, y = u
  α, β, δ, γ = p
  dx = α*x - β*x*y
  dy = -δ*y + γ*x*y
  [dx,dy]
end
function lotka_volterra_noise(du,u,p,t)
  du[1] = 0.01u[1]
  du[2] = 0.01u[2]
end
function lotka_volterra_noise(u,p,t)
  [0.01u[1],0.01u[2]]
end
prob = SDEProblem(lotka_volterra,lotka_volterra_noise,[1.0,1.0],(0.0,10.0))
p = [2.2, 1.0, 2.0, 0.4]
function predict_fd_sde(p)
  solve(prob,SOSRI(),p=p,saveat=0.0:0.1:0.5,sensealg=ForwardDiffSensitivity())[1,:]
end
loss_fd_sde(p) = sum(abs2,x-1 for x in predict_fd_sde(p))
loss_fd_sde(p)

prob = SDEProblem{false}(lotka_volterra,lotka_volterra_noise,[1.0,1.0],(0.0,10.0))
p = [2.2, 1.0, 2.0, 0.4]
function predict_fd_sde(p)
  solve(prob,SOSRI(),p=p,saveat=0.0:0.1:0.5,sensealg=ForwardDiffSensitivity())[1,:]
end
loss_fd_sde(p) = sum(abs2,x-1 for x in predict_fd_sde(p))
loss_fd_sde(p)

@test !iszero(Zygote.gradient(loss_fd_sde,p)[1])

prob = SDEProblem(lotka_volterra,lotka_volterra_noise,[1.0,1.0],(0.0,0.5))
function predict_rd_sde(p)
  solve(prob,SOSRI(),p=p,saveat=0.0:0.1:0.5,sensealg=TrackerAdjoint())[1,:]
end
loss_rd_sde(p) = sum(abs2,x-1 for x in predict_rd_sde(p))
@test !iszero(Zygote.gradient(loss_rd_sde,p)[1])

prob = SDEProblem{false}(lotka_volterra,lotka_volterra_noise,[1.0,1.0],(0.0,0.5))
function predict_rd_sde(p)
  solve(prob,SOSRI(),p=p,saveat=0.0:0.1:0.5,sensealg=TrackerAdjoint())[1,:]
end
loss_rd_sde(p) = sum(abs2,x-1 for x in predict_rd_sde(p))
@test !iszero(Zygote.gradient(loss_rd_sde,p)[1])


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

"""
    LibGit2.GitStatus(repo::GitRepo; status_opts=StatusOptions())

Collect information about the status of each file in the git
repository `repo` (e.g. is the file modified, staged, etc.).
`status_opts` can be used to set various options, for instance
whether or not to look at untracked files or whether to include
submodules or not. See [`StatusOptions`](@ref) for more information.
"""
function GitStatus(repo::GitRepo; status_opts=StatusOptions())
    stat_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_status_list_new, :libgit2), Cint,
                  (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Ptr{StatusOptions}),
                  stat_ptr_ptr, repo.ptr, Ref(status_opts))
    return GitStatus(repo, stat_ptr_ptr[])
end

function Base.length(status::GitStatus)
    return Int(ccall((:git_status_list_entrycount, :libgit2), Csize_t,
                      (Ptr{Ptr{Cvoid}},), status.ptr))
end

function Base.getindex(status::GitStatus, i::Integer)
    1 <= i <= length(status) || throw(BoundsError())
    GC.@preserve status begin
        entry_ptr = ccall((:git_status_byindex, :libgit2),
                          Ptr{StatusEntry},
                          (Ptr{Cvoid}, Csize_t),
                          status.ptr, i-1)
        entry_ptr == C_NULL && throw(Error.GitError(Error.ERROR))
        entry = unsafe_load(entry_ptr)
    end
    return entry
end

"""
    LibGit2.status(repo::GitRepo, path::String) -> Union{Cuint, Cvoid}

Lookup the status of the file at `path` in the git
repository `repo`. For instance, this can be used
to check if the file at `path` has been modified
and needs to be staged and committed.
"""
function status(repo::GitRepo, path::String)
    status_ptr = Ref{Cuint}(0)
    ret =  ccall((:git_status_file, :libgit2), Cint,
                  (Ref{Cuint}, Ptr{Cvoid}, Cstring),
                  status_ptr, repo.ptr, path)
    (ret == Cint(Error.ENOTFOUND) || ret == Cint(Error.EAMBIGUOUS)) && return nothing
    return status_ptr[]
end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
    using Revise, Parameters, Logging, Test
    using CUDA
    includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
    using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
    using ..BitWiseUtils,..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates

    for i in 1:100000
        result = CUDA.zeros(1)
        sourceCPU = rand(UInt32,1)[1]
        targetCPU = rand(UInt32,1)[1]
        source= CuArray([sourceCPU])
        target= CuArray([targetCPU])

        function testKernelPassOnes(source,target,result)
            result[1] = @bitPassOnes(source[1],target[1])
            return
        end
        @cuda threads=32 blocks=1 testKernelPassOnes(source,target,result)
        result[1]
    end


text:
export nn_policy_maker

using Flux: softmax

# an init_transform
function div_by_11(x)
    x ./ 11
end

"""
Given a neural work (nn) that accepts a 4x4x1xn array of inputs and returns
a policy
"""
function nn_policy_maker(nn, init_transform=div_by_11)
    function(board::Bitboard)
        tmp = Float32.(bitboard_to_array(board)) |> init_transform
        reshape(softmax(nn(reshape(tmp, 4, 4, 1, 1))), 4)
    end
end


text:
@testset "NumField/NfAbs" begin
  include("NfAbs/ConjugatesNS.jl")
  include("NfAbs/MPolyGcd.jl")
  include("NfAbs/MPolyAbsFact.jl")
end


text:
using Bokeh
plotfile("jagged_array.html")
hold(true)
# jagged array, this is an array of arrays, not something we should use if possible,
# but Bokeh.jl can deal with it
ys=Vector[[1,2,3], [4,5,6], [7,8,9,10]]
plot(ys, "ro|ko|bo")

x=[14,13,12,11]
plot(x, ys, "rd|kd|bd")

xs = Vector[[7,4,2], [6,5,4], [10.1,10,8, 7.5]]
plot(xs, ys, "rs|ks|bs")
showplot()


text:
@testset "benchmark_normals" begin
    p, q = synthetic_gradient(SynthSphere(50))
    p2, q2 = synthetic_gradient(SynthSphere(51))
    error = benchmark_normals(p, q, p2, q2)
    @test error ≈ 1.4866545112360603
end


text:
##
# Abstract relation.
abstract type AbstractRelation end

##
# Abstract proposotional relation.
abstract type AbstractPropositionalRelation     <: AbstractRelation end
# Abstract modal relation.
abstract type AbstractModalRelation             <: AbstractRelation end

##
# Abstract temporal relation.
abstract type AbstractTemporalRelation          <: AbstractModalRelation end
# Abstract spatial relation.
abstract type AbstractSpatialRelation           <: AbstractModalRelation end
# Abstract spatio-temporal relation.
abstract type AbstractSpatioTemporalRelation    <: AbstractModalRelation end

##
# Abstract point (temporal) relation.
abstract type AbstractPointRelation             <: AbstractTemporalRelation end
# Abstract interval (temporal) relation.
abstract type AbstractIntervalRelation          <: AbstractTemporalRelation end
# Abstract rectangle (spatial) relation.
abstract type AbstractRectangleRelation         <: AbstractSpatialRelation end
# Abstract cube (spatio-temporal) relation.
abstract type AbstractCubeRelation              <: AbstractSpatioTemporalRelation end

##
# Abstract existential interval (temporal) relation.
abstract type AbstractExistentialIntervalRelation   <: AbstractIntervalRelation end
# Abstract universal interval (temporal) relation.
abstract type AbstractUniversalIntervalRelation     <: AbstractIntervalRelation end

##
#
struct IdentityIntervalRelation         <: AbstractIntervalRelation
    algebra::AbstractAlgebra
end
# Using the singleton design pattern; similar to Val{T}.
struct ExistentialIntervalRelation{T}   <: AbstractExistentialIntervalRelation
    algebra::AbstractAlgebra
end
ExistentialIntervalRelation(s::AbstractString,a::AbstractAlgebra)  = ExistentialIntervalRelation{Symbol(s)}(a)
ExistentialIntervalRelation(s::Symbol,a::AbstractAlgebra)          = ExistentialIntervalRelation{s}(a)
struct UniversalIntervalRelation{T}     <: AbstractUniversalIntervalRelation
    algebra::AbstractAlgebra
end
UniversalIntervalRelation(s::AbstractString,a::AbstractAlgebra)    = UniversalIntervalRelation{Symbol(s)}(a)
UniversalIntervalRelation(s::Symbol,a::AbstractAlgebra)            = UniversalIntervalRelation{s}(a)

# TODO disable r.algebra
# TODO extend to higher types (in the hierarchy); e.g., AbstractModalRelation
algebra(r::AbstractIntervalRelation) = getproperty(r, :algebra)

##
Base.show(io::IO, ::IdentityIntervalRelation)           = print(io, "|=|")

Base.show(io::IO, ::ExistentialIntervalRelation{:G})    = print(io, "⟨G⟩")

Base.show(io::IO, ::ExistentialIntervalRelation{:L})    = print(io, "⟨L⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:A})    = print(io, "⟨A⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:O})    = print(io, "⟨O⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:E})    = print(io, "⟨E⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:D})    = print(io, "⟨D⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:B})    = print(io, "⟨B⟩")

Base.show(io::IO, ::ExistentialIntervalRelation{:InvL}) = print(io, "⟨InvL⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:InvA}) = print(io, "⟨InvA⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:InvO}) = print(io, "⟨InvO⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:InvE}) = print(io, "⟨InvE⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:InvD}) = print(io, "⟨InvD⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:InvB}) = print(io, "⟨InvB⟩")

Base.show(io::IO, ::UniversalIntervalRelation{:G})      = print(io, "[G]")

Base.show(io::IO, ::UniversalIntervalRelation{:L})      = print(io, "[L]")
Base.show(io::IO, ::UniversalIntervalRelation{:A})      = print(io, "[A]")
Base.show(io::IO, ::UniversalIntervalRelation{:O})      = print(io, "[O]")
Base.show(io::IO, ::UniversalIntervalRelation{:E})      = print(io, "[E]")
Base.show(io::IO, ::UniversalIntervalRelation{:D})      = print(io, "[D]")
Base.show(io::IO, ::UniversalIntervalRelation{:B})      = print(io, "[B]")

Base.show(io::IO, ::UniversalIntervalRelation{:InvL})   = print(io, "[InvL]")
Base.show(io::IO, ::UniversalIntervalRelation{:InvA})   = print(io, "[InvA]")
Base.show(io::IO, ::UniversalIntervalRelation{:InvO})   = print(io, "[InvO]")
Base.show(io::IO, ::UniversalIntervalRelation{:InvE})   = print(io, "[InvE]")
Base.show(io::IO, ::UniversalIntervalRelation{:InvD})   = print(io, "[InvD]")
Base.show(io::IO, ::UniversalIntervalRelation{:InvB})   = print(io, "[InvB]")

# ##
# # TODO not necessary?
# function check_relation(r::RT, os::NTuple{2,OT}) where {RT,OT}
#     check_relation(OntologicalTrait(RT), OntologicalTrait(OT), r, os)
# end
#
# function check_relation(::IntervalTrait, ::IntervalTrait, r::AbstractIntervalRelation, os::NTuple{2,Interval})
#     if iscrisp(algebra(r))
#         check_crisp_relation(IntervalTrait(), r, os)
#     else
#         check_fuzzy_relation(IntervalTrait(), r, os)
#     end
# end
#
# # TODO it is cumbersome: could be made better?
#
# function check_crisp_relation(::IntervalTrait, r::AbstractIntervalRelation, os::NTuple{2,Interval})
#     check_crisp_relation(r, os)
# end
#
# function check_fuzzy_relation(::IntervalTrait, r::AbstractIntervalRelation, os::NTuple{2,Interval})
#     check_fuzzy_relation(r, os)
# end
#
# # TODO it is cumbersome: os[1].y.x blah..
# # TODO aprire le tuple cosi posso scambiare i due intervalli per le relazioni inverse
#
# function check_crisp_relation(::ExistentialIntervalRelation{:L}, os::NTuple{2,Interval})
#     return os[1].y.x < os[2].x.x ? 1.0 : 0.0
# end
#
# function check_crisp_relation(::ExistentialIntervalRelation{:A}, os::NTuple{2,Interval})
#     return os[1].y.x = os[2].x.x ? 1.0 : 0.0
# end
#
# # function check_crisp_relation(::ExistentialIntervalRelation{:O}, os::NTuple{2,Interval})
# #     return os[1].x.x < os[2].x.x && os[2].x.x < os[1].y.x &&
# # end
#
# function check_fuzzy_relation(::ExistentialIntervalRelation{:L}, os::NTuple{2,Interval})
#     return <̃ₕ(os[1].y.x, os[2].x.x)
# end
#
#
# # my_f(f, a) = "hahaha"
# # my_f(f::typeof(sum), a) = sum(a)
# # my_f(f::typeof(sum), a; gamma=0.7) = "haha !! gamma!!"
# #
# # my_min(f, a) = "default min"
# # my_min(f::typeof(min), a; gamma=0.7) = g(a; gamma)
# #
# # g(a; gamma=0.7) = "min with gamma = $gamma"


text:
"""
    make_sin(n, start, stop; noise = 0.3, f_rand = randn) -> x, y

Generates `n` noisy equally spaces samples of a sinus from `start` to `stop`
by adding `noise .* f_rand(length(x))` to the result of `fun(x)`.
"""
function make_sin(n::Int = 50, start::Real = 0, stop::Real = 2π; noise::Real = 0.3, f_rand::Function = randn)
    x = collect(range(start, stop=stop, length=n))
    y = sin.(x) .+ noise  .* f_rand.()
    return x, y
end

"""
    make_poly(coef, x; noise = 0.01, f_rand = randn) -> x, y

Generates a noisy response for a polynomial of degree `length(coef)`
using the vector `x` as input and adding `noise .* f_randn(length(x))` to the result.
The vector `coef` contains the coefficients for the terms of the polynome.
The first element of `coef` denotes the coefficient for the term with
the highest degree, while the last element of `coef` denotes the intercept.
"""
function make_poly(coef::AbstractVector{R}, x::AbstractVector{T};
                noise::Real = 0.1, f_rand::Function = randn) where {T<:Real,R<:Real}
    n = length(x)
    m = length(coef)
    x_vec = collect(x)
    y = zeros(n)
    @inbounds for i = 1:n
        for k = 1:m
            y[i] += coef[k] * x_vec[i]^(m-k)
        end
    end
    y .+= noise .* f_rand.()
    x_vec, y
end

"""
    make_spiral(n, a, theta, b; noise = 0.01, f_rand = randn) -> x, y

Generates `n` noisy responses for a spiral with two labels. Uses the radius, angle
and scaling arguments to space the points in 2D space and adding `noise .* f_randn(n)`
to the response.
"""
function make_spiral(n::Int = 97, a::Real = 6.5, theta::Real = 16.0, b::Real=104.0;
                noise::Real = 0.1, f_rand::Function = randn)
    x = zeros(Float64, (2, 2*n))
    y = zeros(Int, 2*n)
    index = 0:1.0:(n-1)
    for i = 1:n
        _angle = index[i]*pi/theta
    	_radius = a * (b-index[i]) / b
    	x_coord = _radius * sin(_angle)
    	y_coord = _radius * cos(_angle)
    	x[1, i] = x_coord
    	x[2, i] = y_coord
        x[1, n+i] = -(x_coord)
    	x[2, n+i] = -(y_coord)
        y[i] = 1
        y[n+i] = 0
    end
    x[1, :] .+= noise .* f_rand.()
    x[2, :] .+= noise .* f_rand.()
    x, y
end


text:
using DataFrames, CSV

function mwe()
  df = DataFrame(rand(1000,5))

  #seems to work without the csv usage

  println(Threads.nthreads())
  dfnocsv = deepcopy(df)
  filter!(r->r.x1<0.5, dfnocsv)
  println("size dfnocsv: ", size(dfnocsv))

  #sending the df to a CSV first breaks the code
  df |> CSV.write("test.csv")
  dfcsv = CSV.read("test.csv", copycols=false, threaded=false) |> DataFrame
  println(typeof(dfcsv.x1))
  filter!(r->r.x1<0.5, dfcsv)
  println("size dfcsv: ", size(dfcsv))

end

mwe()

@inline function fixcsv(df::DataFrame)::DataFrame
  #=for c::Int ∈ 1:size(df,2)
    df[!,c] = Array(df[!,c])
  end=#
  return DataFrame(view(df, :, :))
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Dates
using Random

include("init.jl")
include("dynamics.jl")
include("analysis.jl")

#Runs the simulation at a determined temperature for a determined number of steps
function Run(temperature,nodes,steps,P)
    global T=temperature
    global J = 1

    Random.seed!(1234)

    #Get time of simulation
    exportTime = Dates.format(Dates.now(), "yyyy-mm-ddTHH-MM-SS")

    #Initializing network and data frame
    createDataFrame()
    initNetwork(nodes)

    #Executing steps
    iter = div(steps,100)
    for i in 1:iter
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)
        Procedure2(rand(1:nodes),nodes,P)

        if mod(i*100, div(steps, 100)) == 0
            println(div(i*10000,steps),"% done (",i*100,")")
        end
    end

    #Export data to a new folder
    exportData(exportTime)
    plotAnalysis(steps, exportTime)
end

#Runs the simulation at T for a determined time interval (in hours)
function SetRun(temperature,nodes,duration)
    global T=temperature
    global J = 1

    #Get time of simulation
    exportTime = Dates.format(Dates.now(), "yyyy-mm-ddTHH-MM-SS")

    #Initializing network and data frame
    createDataFrame()
    initNetwork(nodes)

    timelimit = Dates.now() + Dates.Minute(duration)

    #Executing steps
    while Dates.now() < timelimit #&& length(Data.E) <= 2000000
        Procedure2(rand(1:nodes),nodes,0.5)
    end

    #Export data to a new folder
    exportData(exportTime)
    plotAnalysis(length(Data.E)-1, exportTime)
end

#Runs n determined length simulation with populations increasing 100 agents at a time up to a determined size
function SweepPopulations(temperature,maxNodes,steps)
    global T=temperature
    global J = 1

    #Get time of simulation
    #exportTime = Dates.format(Dates.now(), "yyyy-mm-ddTHH-MM-SS")

    for i in 100:100:maxNodes
        #Initializing network and data frame
        createDataFrame()
        initNetwork(i)

        #Executing steps
        for j in 1:steps
            Procedure2(rand(1:nodes),nodes,P)
        end

        #Export data to a new folder
        exportData(string("Population Sweep: ",i))
        plotAnalysis(steps, string("Population Sweep: ",i))
    end
end

#Run successive simulations in increments of 100 steps
function SweepTime(temperature,nodes,maxSteps)
    global T=temperature
    global J = 1

    #Get time of simulation
    #exportTime = Dates.format(Dates.now(), "yyyy-mm-ddTHH-MM-SS")

    for i in 100:100:maxSteps
        #Initializing network and data frame
        createDataFrame()
        initNetwork(nodes)

        #Executing steps
        for j in 1:i
            Procedure2(rand(1:nodes),nodes,0.5)
        end

        #Export data to a new folder
        exportData(string("Time Sweep: ",i))
        plotAnalysis(i, string("Time Sweep: ",i))
    end
end

#Runs multiple simulation at incremental temperatures
function SweepTemperature(T_min,T_max,nodes,steps,P)

end

max_stars_repo_path: Scripts/Prototype 2/main.jl

max_stars_repo_name: luca-mm/Bachelor-Thesis

max_issues_repo_path: Scripts/Prototype 2/main.jl

max_issues_repo_name: luca-mm/Bachelor-Thesis

max_forks_repo_path: Scripts/Prototype 2/main.jl

max_forks_repo_name: luca-mm/Bachelor-Thesis


text:
function flip(point::AbstractArray, prob::Float64)
    @assert zero(prob)<=prob<=one(prob) "The probability must be in range [0,1]"
    α = deepcopy(point)
    if prob == zero(prob)
        return α
    elseif prob == one(prob)
        return  flip.(α)
    end
    dist = Geometric(prob)
    index = rand(dist)+1
    while index <= length(α)
        α[index] = flip(α[index])
        index = index + rand(dist)+1
    end
    α
end

flip(bit::Bool) = ~bit

function staticBitMutation(prob::Float64)
    function _Mutation(instance, candidate, fitness)
        tmp = flip(candidate.point, prob)
        (point=tmp, fitness=fitness(instance, tmp))
    end
    _Mutation
end

max_stars_repo_name: ahrajabi/DiscO

max_issues_repo_name: ahrajabi/DiscO

max_forks_repo_name: ahrajabi/DiscO


text:
@discord_object struct Presence
    user::User
    roles::Vector{Snowflake}
    game::Activity
    guild_id::Snowflake
    status::String
    activities::Vector{Activity}
    client_status::ClientStatus
    premium_since::Union{String, DateTime}
    nick::String
end


text:
function include_all_files(cur_item; package_name=nothing, is_testing=false, is_focused=false, is_sorted=false, reload_function=Nullable, is_revised=true)
  loaded_files = []
  all_files = get_all_files(cur_item, package_name=package_name, is_testing=is_testing, is_sorted=is_sorted)
  unloaded_files = setdiff(all_files, loaded_files)

  cur_module = getfield(Main, Symbol(package_name))

  while length(unloaded_files) > 0
    new_file_count = 0

    for file in all_files
      is_already_loaded = in(file, loaded_files)
      if is_already_loaded ; continue ; end

      init_methods_list = filter(
        cur_list -> !isempty(cur_list),
        map(
          cur_var -> try methods(getfield(cur_module, cur_var)) ; catch [] ; end,
          Julz.get_all_symbols(cur_module)
        )
      )

      try

        if is_testing

          ignore_file = is_focused && !check_for_focus(file)
          ignore_file |= check_for_skip(file)

          if ignore_file
            push!(loaded_files, file)
            new_file_count += 1
            continue
          end

          if reload_function == Nullable
            error("Need to pass reload_function for tests")
          end

          defaults_file_name = "defaults.jl"
          reload_function(defaults_file_name, true)

          input_file_name = "test/input.jl"
          reload_function(input_file_name, true)

        end

        include(file)

        skip_revise = is_testing || !is_revised
        skip_revise |= contains(file, "config/initializers")
        skip_revise |= contains(file, "src/modules")

        if !skip_revise
          Revise.track(cur_module, file)
        end

      catch

        cur_methods_list = filter(
          cur_list -> !isempty(cur_list),
          map(
            cur_var -> try methods(getfield(cur_module, cur_var)) ; catch [] ; end,
            Julz.get_all_symbols(cur_module)
          )
        )

        cur_methods_list = setdiff(cur_methods_list, init_methods_list)

        cur_methods_list = map(
          cur_methods -> filter(
            cur_method -> cur_method.module == cur_module,
            cur_methods.ms
          ),
          cur_methods_list
        )

        cur_methods = flatten(cur_methods_list)

        for cur_method in cur_methods
          delete_method(cur_method)
        end

        Docs.initmeta(cur_module)

        for cur_method in cur_methods
          cur_method_name = cur_method.name

          b = Docs.Binding(cur_module, cur_method_name)
          m = get!(Docs.meta(cur_module), b, Docs.MultiDoc())

          cur_sig_string = string(cur_method.sig)

          cur_sig_string = replace(cur_sig_string, "$(cur_module).#$(cur_method_name)", "")
          cur_sig_string = replace(cur_sig_string, "{,", "{")

          cur_sig = eval(parse(string(cur_sig_string)))

          cur_keys = collect(keys(m.docs))

          for cur_key in cur_keys
            ( cur_sig <: cur_key ) || continue

            delete!(m.docs, cur_key)
          end
        end

        continue
      end

      push!(loaded_files, file)
      new_file_count += 1
    end

    unloaded_files = setdiff(all_files, loaded_files)

    if new_file_count == 0
      bad_file = unloaded_files[1]
      println(bad_file)
      include(bad_file)
    end
  end
end


text:
module TestPatterns
using Test

using ReTest: and, or, not, interpolated

@testset "patterns: ==" begin
    basics = [and(), or(), not(0), interpolated, 0, r""]
    for a = basics, b = basics
        if a === b
            @test a == b
        else
            @test a != b
            for f in (and, or)
                @test f(a) == f(a)
                @test f(a) != f(b)
                @test f(a, b) == f(a, b)
                @test f(a, b) != f(b, a)
                @test f(and(a, b), or(b, a)) == f(and(a, b), or(b, a))
                @test f(and(a, b), or(a, b)) != f(and(a, b), or(b, a))
                @test f(not(a), and(a, b)) == f(not(a), and(a, b))
                @test f(and(a, b), not(a)) != f(not(a), and(a, b))
                @test f(not(a), r"a", and(a, b)) == f(not(a), r"a", and(a, b))
                @test f(not(a), r"a", and(a, b)) != f(not(a), r"b", and(a, b))
            end
            @test not(a) == not(a)
            @test not(a) != not(b)
            @test not(not(a)) == not(not(a))
        end
    end
end

end # TestPatterns


text:
"""
Jump dispersal rules.
"""
abstract type AbstractJumpDispersal{R,W} <: ManualRule{R,W} end

"""
    JumpDispersal(spotrange)
    JumpDispersal{R,W}(spotrange)
    JumpDispersal(; spotrange=30.0)

Jump dispersal simulates random long distance dispersal events. A random cell within
the `spotrange` is invaded.

Pass grid name `Symbol`s to `R` and `W` type parameters to use specific grids.

$(FIELDDOCTABLE)
"""
@Probabilistic struct JumpDispersal{R,W,SR} <: AbstractJumpDispersal{R,W}
    # Field       | Def  | Flatten | Limits       | Description
    spotrange::SR | 30.0 | true    | (0.0, 100.0) | "A number or Unitful.jl distance with the same units as cellsize"
end

# TODO update this and test
@inline applyrule!(data, rule::AbstractJumpDispersal{R,W}, state, index
                  ) where {R,W} = begin
    # Ignore empty cells
    state > zero(state) || return state

    # Random dispersal events
    rand() < rule.prob_threshold || return state

    # Randomly select spotting distance
    intspot = round(Int, rule.spotrange)
    rnge = -intspot:intspot
    jump = (rand(rnge), rand(rnge))
    jumpdest, is_inbounds = inbounds(jump .+ index, gridsize(data), RemoveOverflow())

    # Update spotted cell if it's on the grid
    if is_inbounds
        @inbounds data[W][jumpdest...] = state
    end

    state
end


text:
"""
    opening(img; dims=coords_spatial(img), r=1)
    opening(img, se)

Perform the morphological opening on `img`. The opening operation is defined as erosion
followed by a dilation: `dilate(erode(img, se), se)`.

$(_docstring_se)

# Examples

```jldoctest; setup = :(using ImageMorphology)
julia> img = trues(7,7); img[2, 2] = false; img[3:5, 3:5] .= false; img[4, 4] = true; img
7×7 BitMatrix:
 1  1  1  1  1  1  1
 1  0  1  1  1  1  1
 1  1  0  0  0  1  1
 1  1  0  1  0  1  1
 1  1  0  0  0  1  1
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1

julia> opening(img)
7×7 BitMatrix:
 0  0  1  1  1  1  1
 0  0  1  1  1  1  1
 1  1  0  0  0  1  1
 1  1  0  0  0  1  1
 1  1  0  0  0  1  1
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1

julia> opening(img, strel_diamond(img)) # use diamond shape SE
7×7 BitMatrix:
 1  1  1  1  1  1  1
 1  0  1  1  1  1  1
 1  1  0  0  0  1  1
 1  1  0  0  0  1  1
 1  1  0  0  0  1  1
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1
```

## See also

- [`opening!`](@ref) is the in-place version of this function.
- [`closing`](@ref) is the dual operator of `opening` in the sense that
  `complement.(opening(img)) == closing(complement.(img))`.
"""
opening(img; kwargs...) = opening!(similar(img), img, similar(img); kwargs...)
opening(img, se) = opening!(similar(img), img, se, similar(img))

"""
    opening!(out, img, buffer; [dims], [r])
    opening!(out, img, se, buffer)

The in-place version of [`opening`](@ref) with input image `img` and output image `out`. The
intermediate erosion result is stored in `buffer`.
"""
function opening!(out, img, buffer; dims=coords_spatial(img), r=nothing)
    return opening!(out, img, strel_box(img, dims; r), buffer)
end
function opening!(out, img, se, buffer)
    erode!(buffer, img, se)
    dilate!(out, buffer, se)
    return out
end
opening!(out::AbstractArray{<:Color3}, img, se, buffer) = throw(ArgumentError("color image is not supported"))


text:
module TurnsLogger

using ...PCG.Topologies: Topology
using ...PCG.Types: Recorder
using ...PCG.Storages
using ...PCG.Universes
using ...PCG.Universes: Universe

export TurnsLoggerRecorder

const NO_TURNS_LIMIT = -1::Int64

struct TurnsLoggerRecorder <: Recorder
    total_turns::Int64

    TurnsLoggerRecorder(total_turns::Int64=NO_TURNS_LIMIT) = new(total_turns)
end

function Universes.record_state!(universe::Universe, recorder::TurnsLoggerRecorder)
    if recorder.total_turns == NO_TURNS_LIMIT
        println("turn $(universe.turn) processed")
    else
        println("turn $(universe.turn)/$(recorder.total_turns) processed")
    end
end

function Universes.finish_recording!(recorder::TurnsLoggerRecorder)
    println("processed")
end

end

max_stars_repo_name: Tiendil/pcg-julia

max_issues_repo_name: Tiendil/pcg-julia

max_forks_repo_name: Tiendil/pcg-julia


text:
using Base:GenericIOBuffer
import Base:read, write

include("./typemaps.jl")
include("./utils.jl")

const SRID_FLAG = UInt32(0x20)

@enum Endian BigEndian=0x00 LittleEndian=0x01

function sameendinan(flag::UInt8)
    if ENDIAN_BOM == 0x01020304 && Endian(flag) == BigEndian
        true
    elseif ENDIAN_BOM == 0x04030201 && Endian(flag) == LittleEndian
        true
    else
        false
    end
end

function sameendinan(endian::Endian)
    if ENDIAN_BOM == 0x01020304 && endian == BigEndian
        true
    elseif ENDIAN_BOM == 0x04030201 && endian == LittleEndian
        true
    else
        false
    end
end

"""
    WKB Type
"""
struct WKB
    data::AbstractVector{UInt8}
end

function WKB(s::AbstractString; hex=false)
    hex ? WKB(hex2bytes(s)) : WKB(Vector{UInt8}(s))
end

macro WKB_str(s)
    WKB(s)
end

macro HEXWKB_str(s)
    WKB(s; hex=true)
end

IOBuffer(wkb::WKB) = Base.IOBuffer(wkb.data)

read(s::IO, T; swap=false) = swap ? read(s, T) |> bswap : read(s, T)
write(s::IO, T; swap=false) = swap ? write(s, T |> bswap) : write(s, T)

"""
    read(WKB)
"""
function read(b::WKB; withsrid=true, zipmulti=true)
    wkbio = IOBuffer(b)
    _read(wkbio, withsrid=withsrid, zipmulti=zipmulti)
end

function write(geo::NamedTuple; endian=LittleEndian, hex=false)
    fields = keys(geo)
    if :type in fields
        geotype = geo.type
    else
        error("Geometry object must have type fields?")
    end

    if :coordinates in fields
        coords_or_geoms = geo.coordinates
    elseif :geometries in fields
        coords_or_geoms = geo.geometries
    else
        error("Geometry object must have coordinates or geometries fields?")
    end

    meta = :meta in fields ? geo.meta : nothing
    srid = meta != nothing && :srid in keys(meta) ? meta.srid : nothing

    io = Base.IOBuffer()
    io = _write(io, endian=endian)
    bytes = take!(io)
    hex ? String(bytes2hex(bytes)) : String(bytes)
end

function _write(io::GenericIOBuffer, endian::Endian, typecode::UInt32, srid::UInt32=nothing)
    write(io, endian |> UInt8)
    swap = !sameendinan(endian)

    if srid != nothing
        typefield = (typecode) | bswap(SRID_FLAG)
        write(typecode, swap=swap)
        write(srid, swap=swap)
    else
        typefield = typecode
        write(typecode, swap=swap)
    end
end

function _read(io::GenericIOBuffer; withsrid=true, onlydata=false, zipmulti=true)
    endian = read(io, UInt8)
    srid = nothing
    same = sameendinan(endian)

    typefield = read(io, UInt32, swap=!same)
    typecode = typefield & 0x00ffffff
    geotype = BINARY_TO_GEOM_TYPE[typecode] |> string
    hassrid = (typefield & 0xff000000) |> bswap == SRID_FLAG
    if hassrid
        srid = read(io, UInt32, swap=!same)
    end

    if startswith(geotype, "Multi")
        coords = _load(WKBCode(typecode), io, swap=!same, onlydata=zipmulti)
    else
        coords = _load(WKBCode(typecode), io, swap=!same)
    end

    if onlydata
        coords
    elseif hassrid && withsrid
        if geotype == "GeometryCollection"
            (type=geotype, geometries=coords,
             meta=(srid=srid,),
             crs=(type="name", properties=(name="EPSG$srid",)))
        else
            (type=geotype, coordinates=coords,
             meta=(srid=srid,),
             crs=(type="name", properties=(name="EPSG$srid",)))
        end
    else
        if geotype == "GeometryCollection"
            (type=geotype, geometries=coords)
        else
            (type=geotype, coordinates=coords)
        end
    end
end

"""
    readcoords(data::GenericIOBuffer, n=2, swap=false)
    n坐标的维度
"""
function readcoords(io::GenericIOBuffer; nb=2, swap=false)
    [read(io, Float64, swap=swap) for i in 1:nb]
end

function writecoords(io::GenericIOBuffer, data::Vector{Float64}; swap=false)
    nb = length(data)
    for i in 1:nb
        write(io, data[i], swap=swap)
    end
end

function readnodes(io::GenericIOBuffer; nb=2, swap=false)
    node = read(io, UInt32, swap=swap)
    [readcoords(io, nb=nb, swap=swap) for i in 1:node]
end

function writenodes(io::GenericIOBuffer, data::Vector{Vector{Float64}}; swap=false)
    node = length(data) |> UInt32
    write(io, node, swap=swap)
    for i in 1:node
        writecoords(io, data[i], swap=swap)
    end
end

function readrings(io::GenericIOBuffer; nb=2, swap=false)
    ring = read(io, UInt32, swap=swap)
    [readnodes(io, nb=nb, swap=swap) for i in 1:ring]
end

function writerings(io::GenericIOBuffer, data::Vector{Vector{Vector{Float64}}}; swap=false)
    ring = length(data) |> UInt32
    write(io, ring, swap=swap)
    for i in 1:ring
        writenodes(io, data[i], swap=swap)
    end
end

function readmulti(io::GenericIOBuffer; swap=false, onlydata=true)
    num = read(io, UInt32, swap=swap)
    [_read(io; onlydata=onlydata) for i in 1:num]
end

function writemulti(io::GenericIOBuffer, data::Vector{Vector{Vector{Vector{Float64}}}}; swap=false)
    num = length(data) |> UInt32
    write(io, num, swap=swap)

end

function readgc(io::GenericIOBuffer; swap=false, onlydata=false)
    num = read(io, UInt32, swap=swap)
    [_read(io; onlydata=onlydata) for i in 1:num]
end

"""
根据WKB type code值派发
"""
struct WKBCode{x}
end

WKBCode(x) = (Base.@_pure_meta; WKBCode{x}())

function _load(::WKBCode{WKB_2D[:Point]}, data::GenericIOBuffer; swap=false)
    readcoords(data, swap=swap)
end

function _load(::WKBCode{WKB_Z[:Point]}, data::GenericIOBuffer; swap=false)
    readcoords(data, nb=3, swap=swap)
end

function _load(::WKBCode{WKB_M[:Point]}, data::GenericIOBuffer; swap=false)
    readcoords(data, nb=3, swap=swap)
end

function _load(::WKBCode{WKB_ZM[:Point]}, data::GenericIOBuffer; swap=false)
    readcoords(data, nb=4, swap=swap)
end

function _load(::WKBCode{WKB_2D[:LineString]}, data::GenericIOBuffer; swap=false)
    readnodes(data, swap=swap)
end

function _load(::WKBCode{WKB_Z[:LineString]}, data::GenericIOBuffer; swap=false)
    readnodes(data, nb=3, swap=swap)
end

function _load(::WKBCode{WKB_M[:LineString]}, data::GenericIOBuffer; swap=false)
    readnodes(data, nb=3, swap=swap)
end

function _load(::WKBCode{WKB_ZM[:LineString]}, data::GenericIOBuffer; swap=false)
    readnodes(data, nb=4, swap=swap)
end

function _load(::WKBCode{WKB_2D[:Polygon]}, data::GenericIOBuffer; swap=false)
    readrings(data, swap=swap)
end

function _load(::WKBCode{WKB_Z[:Polygon]}, data::GenericIOBuffer; swap=false)
    readrings(data, nb=3, swap=swap)
end

function _load(::WKBCode{WKB_M[:Polygon]}, data::GenericIOBuffer; swap=false)
    readrings(data, nb=3, swap=swap)
end

function _load(::WKBCode{WKB_ZM[:Polygon]}, data::GenericIOBuffer; swap=false)
    readrings(data, nb=4, swap=swap)
end

function _load(::WKBCode{WKB_2D[:MultiPoint]}, data::GenericIOBuffer; swap=false, onlydata=true)
    readmulti(data; swap=swap, onlydata=onlydata)
end

function _load(::WKBCode{WKB_Z[:MultiPoint]}, data::GenericIOBuffer; swap=false, onlydata=true)
    readmulti(data; swap=swap, onlydata=onlydata)
end

function _load(::WKBCode{WKB_M[:MultiPoint]}, data::GenericIOBuffer; swap=false, onlydata=true)
    readmulti(data; swap=swap, onlydata=onlydata)
end

function _load(::WKBCode{WKB_ZM[:MultiPoint]}, data::GenericIOBuffer; swap=false, onlydata=true)
    readmulti(data; swap=swap, onlydata=onlydata)
end

function _load(::WKBCode{WKB_2D[:MultiLineString]}, data::GenericIOBuffer; swap=false, onlydata=true)
    readmulti(data; swap=swap, onlydata=onlydata)
end

function _load(::WKBCode{WKB_Z[:MultiLineString]}, data::GenericIOBuffer; swap=false, onlydata=true)
    readmulti(data; swap=swap, onlydata=onlydata)
end

function _load(::WKBCode{WKB_M[:MultiLineString]}, data::GenericIOBuffer; swap=false, onlydata=true)
    readmulti(data; swap=swap, onlydata=onlydata)
end

function _load(::WKBCode{WKB_ZM[:MultiLineString]}, data::GenericIOBuffer; swap=false, onlydata=true)
    readmulti(data; swap=swap, onlydata=onlydata)
end

function _load(::WKBCode{WKB_2D[:MultiPolygon]}, data::GenericIOBuffer; swap=false, onlydata=true)
    readmulti(data; swap=swap, onlydata=onlydata)
end

function _load(::WKBCode{WKB_Z[:MultiPolygon]}, data::GenericIOBuffer; swap=false, onlydata=true)
    readmulti(data; swap=swap, onlydata=onlydata)
end

function _load(::WKBCode{WKB_M[:MultiPolygon]}, data::GenericIOBuffer; swap=false, onlydata=true)
    readmulti(data; swap=swap, onlydata=onlydata)
end

function _load(::WKBCode{WKB_ZM[:MultiPolygon]}, data::GenericIOBuffer; swap=false, onlydata=true)
    readmulti(data; swap=swap, onlydata=onlydata)
end

function _load(::WKBCode{WKB_2D[:GeometryCollection]}, data::GenericIOBuffer; swap=false, onlydata=false)
    readgc(data; swap=swap, onlydata=onlydata)
end

function _load(::WKBCode{WKB_Z[:GeometryCollection]}, data::GenericIOBuffer; swap=false, onlydata=false)
    readgc(data; swap=swap, onlydata=onlydata)
end

function _load(::WKBCode{WKB_M[:GeometryCollection]}, data::GenericIOBuffer; swap=false, onlydata=false)
    readgc(data; swap=swap, onlydata=onlydata)
end

function _load(::WKBCode{WKB_ZM[:GeometryCollection]}, data::GenericIOBuffer; swap=false, onlydata=false)
    readgc(data; swap=swap, onlydata=onlydata)
end


text:
using MarketData

facts("time series methods") do

  context("lag takes previous day and timestamps it to next day") do
      @fact lag(cl).values[1]    => roughly(111.94)
      @fact lag(cl).timestamp[1] => Date(2000,1,4)
  end

  context("lag accepts kwarg") do
      @fact lag(cl, period=9).timestamp[1] => Date(2000,1,14)
  end

  context("lag operates on 2d arrays") do
      @fact lag(ohlc, period=9).timestamp[1] => Date(2000,1,14)
  end

  context("lag returns 1d from 1d time arrays") do
      @fact ndims(lag(cl).values) => 1
  end

  context("lag returns 2d from 2d time arrays") do
      @fact ndims(lag(ohlc).values) => 2
  end

  context("lead takes next day and timestamps it to current day") do
      @fact lead(cl).values[1]    => roughly(102.5)
      @fact lead(cl).timestamp[1] => Date(2000,1,3)
  end

  context("lead accepts kwarg") do
      @fact lead(cl, period=9).values[1]    => 100.44
      @fact lead(cl, period=9).timestamp[1] => Date(2000,1,3)
  end

  context("lead operates on 2d arrays") do
      @fact lead(ohlc, period=9).timestamp[1] => Date(2000,1,3)
  end

  context("lead returns 1d from 1d time arrays") do
      @fact ndims(lead(cl).values) => 1
  end

  context("lead returns 2d from 2d time arrays") do
      @fact ndims(lead(ohlc).values) => 2
  end

  context("correct simple return value") do
      @fact percentchange(cl).values[1] => roughly((102.5-111.94)/111.94)
  end

  context("correct log return value") do
      @fact percentchange(cl, method="log").values[1] => roughly(log(102.5) - log(111.94))
  end

  context("moving supplies correct window length") do
      @fact moving(cl, mean, 10).values[1]    => roughly(sum(cl.values[1:10])/10)
      @fact moving(cl, mean, 10).timestamp[1] => Date(2000,1,14)
  end

  context("upto method accumulates") do
      @fact upto(cl, sum).values[10]    => roughly(sum(cl.values[1:10]))
      @fact upto(cl, mean).values[10]   => roughly(sum(cl.values[1:10])/10)
      @fact upto(cl, sum).timestamp[10] => Date(2000,1,14)
  end
end

facts("base element-wise operators on TimeArray values") do

  context("correct alignment and operation between two TimeVectors") do
     @fact (cl .+ op).values[1]           => roughly(216.82)
     @fact (cl .- op).values[1]           => roughly(7.06)
     @fact (cl .* op).values[1]           => roughly(11740.27)
     @fact (cl ./ op).values[1]           => roughly(1.067315)
  end

  context("only values on intersecting Dates computed") do
     @fact (cl[1:2] ./ op[2:3]).values[1] => roughly(0.946882)
     @fact (cl[1:4] .+ op[4:7]).values[1] => roughly(201.12)
     @fact length(cl[1:2] ./ op[2:3])     => 1
     @fact length(cl[1:4] .+ op[4:7])     => 1
  end

  context("correct dot operation between TimeVectors values and Int/Float64 and viceversa") do
     @fact (cl .- 100).values[1] => roughly(11.94)
     @fact (cl .+ 100).values[1] => roughly(211.94)
     @fact (cl .* 100).values[1] => roughly(11194)
     @fact (cl ./ 100).values[1] => roughly(1.1194)
     @fact (cl .^ 2).values[1]   => roughly(12530.5636)
     @fact (100 .- cl).values[1] => roughly(-11.94)
     @fact (100 .+ cl).values[1] => roughly(211.94)
     @fact (100 .* cl).values[1] => roughly(11194)
     @fact (100 ./ cl).values[1] => roughly(0.8933357155619082)
     @fact (2 .^ cl).values[1]   => 4980784073277740581384811358191616
  end

  # these methods need to be deprecated to match base
  context("throw error when calling  non-dot operation between TimeVectors values and Int/Float64 and viceversa") do
     @fact_throws (cl - 100).values[1]
     @fact_throws (cl + 100).values[1]
     @fact_throws (cl * 100).values[1]
     @fact_throws (cl / 100).values[1]
     @fact_throws (cl ^ 2).values[1]
     @fact_throws (100 - cl).values[1]
     @fact_throws (100 + cl).values[1]
     @fact_throws (100 * cl).values[1]
     @fact_throws (100 / cl).values[1]
     @fact_throws (2 ^ cl).values[1]
  end

  context("correct operation between two TimeVectors values returns bool for comparisons") do
     @fact (cl .> op).values[1]  => true
     @fact (cl .< op).values[1]  => false
     @fact (cl .<= op).values[1] => false
     @fact (cl .>= op).values[1] => true
     @fact (cl .== op).values[1] => false
  end

  context("correct operation between TimeVectors values and Int/Float64 (and viceversa) returns bool for comparison") do
     @fact (cl .> 111.94).values[1]  => false
     @fact (cl .< 111.94).values[1]  => false
     @fact (cl .>= 111.94).values[1] => true
     @fact (cl .<= 111.94).values[1] => true
     @fact (cl .== 111.94).values[1] => true
     @fact (111.94 .> cl).values[1]  => false
     @fact (111.94 .< cl).values[1]  => false
     @fact (111.94 .>= cl).values[1] => true
     @fact (111.94 .<= cl).values[1] => true
     @fact (111.94 .== cl).values[1] => true
  end
end

facts("basecall works with Base methods") do

  context("cumsum works") do
    @fact basecall(cl, cumsum).values[2] => cl.values[1] + cl.values[2]
  end

  context("log works") do
    @fact basecall(cl, log).values[2] => log(cl.values[2])
  end
end


text:
mutable struct ConstGate <: QGate
	kernel::GateKernel
	ConstGate(kernel::GateKernel) = new(kernel)
end #struct

function ConstGate(data::Array{T} where {T<:Number}, name::String="Anonymous")
	shape = size(data)
	for i ∈ shape
		if i != 2
			throw(ArgumentError("Given data array is not a valid gate."))
		end
	end
	f() = data
	kernel = GateKernel(f, name)
	return ConstGate(kernel)
end

ConstGate(f::Function, name::String="Anonymous") = ConstGate(GateKernel(f, name))

name(gate::ConstGate) = name(gate.kernel)
data(gate::ConstGate) = func(gate)()
func(gate::ConstGate) = func(gate.kernel)

copy(gate::ConstGate) = gate
==(gate1::ConstGate, gate2::ConstGate) = (gate1.kernel == gate2.kernel)
control(gate::ConstGate) = ConstGate(control(gate.kernel))

function show(io::IO, gate::ConstGate)
	printstyled(io, name(gate); bold=true, color= :blue)
end

inverse(gate::ConstGate) = ConstGate(inverse(gate.kernel))
function randomConstGate(n::Int=rand(1:3))
	if n==1
		choices = [X,Y,Z,H,T,TDAG,S,SDAG]
	elseif n==2
		choices = [CNOT, SWAP]
	else
		choices = [TOFFOLI, FREDKIN]
	end
	count = length(choices)
	return choices[rand(1:count)]
end


text:
using VTUFileHandler
using Test

@testset "VTUFileHandler test" begin
    file = joinpath(".","pointheatsource.vtu")
    tempfieldstr = "temperature_interpolated"
    set_interpolation_keywords([tempfieldstr])
	vtu = VTUFileHandler.VTUFile(file);
	tmp = zero(vtu)
	tmp += vtu
	vtu += tmp
    @test all(vtu[tempfieldstr] .== 2.0*tmp[tempfieldstr])
    vtu -= tmp
    @test all(vtu[tempfieldstr] .== tmp[tempfieldstr])
end

max_stars_repo_name: baxmittens/VTUFileHandler

max_issues_repo_name: baxmittens/VTUFileHandler

max_forks_repo_name: baxmittens/VTUFileHandler


text:
@testset "Constructor tests" begin
    @testset "Empty constructor" begin
        x = SummarizedExperiment()

#        @test x.nrow == 0
#        @test x.ncol == 0
        @test length(x.assays) == 0
        @test size(x.rowdata) == (0, 1)
        @test size(x.coldata) == (0, 1)
        @test length(x.metadata) == 0
    end

    assays = OrderedDict{String, AbstractArray}(
       "foobar" => [[1,2] [3,4] [5,6]],
       "whee" => [[1.2,2.3] [3.4,4.5] [5.6,7.8]])

    @testset "Default arguments" begin
        x = SummarizedExperiment(assays)
#        @test x.nrow == 2
#        @test x.ncol == 3
        @test length(x.assays) == 2
        @test size(x.rowdata) == (2, 1)
        @test size(x.coldata) == (3, 1)
        @test length(x.metadata) == 0
    end

    rowdata = DataFrame(name = ["Gene1", "Gene2"],
                        Features = ["Feat1", "Feat2"])
    coldata = DataFrame(name = ["Sample1", "Sample2", "Sample3"],
                        Treatment = ["normal", "drug1", "drug2"])
    metadata = Dict{String,Any}("foo" => "bar")

    @testset "Full arguments" begin
        x = SummarizedExperiment(assays, rowdata, coldata, metadata)
#        @test x.nrow == 2
#        @test x.ncol == 3
        @test length(x.assays) == 2
        @test size(x.rowdata) == (2, 2)
        @test size(x.coldata) == (3, 2)
        @test length(x.metadata) == 1
    end
end