text:
module ELabFTW

using Dates
using HTTP
using JSON3
using StructTypes

function __init__()
  endpoint = try
    ENV["ELABFTW_ENDPOINT"]
  catch
    ""
  end

  token = try
    ENV["ELABFTW_API_KEY"]
  catch
    ""
  end

  global BASE_API = BaseAPI(
		endpoint = endpoint,
	 	token = token
	)
end

include("request.jl")
include("experiment.jl")
include("item.jl")
include("template.jl")
include("api.jl")

export Experiment, Item, Template
export create_experiment, create_item, create_template
export get_all_experiments, get_experiment
export get_all_items, get_item, get_item_types
export get_all_templates, get_template
export get_tags, get_upload, get_status
export post_experiment, post_item, post_template
export add_link_to_experiment, add_link_to_item
export upload_to_experiment, upload_to_item
export add_tag_to_experiment, add_tag_to_item

end #module


text:
export FORCES,
       ENERGIES,
       VIRIALS,
       nonbonded_computation_tiles,
       compute_nonbonded!,
       naively_compute_nonbonded!

import CUDA

CUDA.allowscalar(false)

const FORCES = 1 << 0
const ENERGIES = 1 << 1
const VIRIALS = 1 << 2

const WARPSIZE = 32

function nonbonded_computation_tiles(N)
    n = Int32(cld(N, WARPSIZE))
    pairs = Vector{Tuple{Int32,Int32}}(undef, n*(n+1)÷2)
    k = 0
    for i=0:n-1, j=1:n-i
        pairs[k+=1] = (j, j+i)
    end
    return CUDA.cu(pairs)
end

function size_adjusted(array::Vector{T}) where {T}
    lacking = WARPSIZE*cld(length(array), WARPSIZE) - length(array)
    return vcat(array, Vector{T}(undef, lacking))
end

function size_adjusted(array::Matrix{T}, transpose=false) where {T}
    rows, cols = transpose ? (2, 1) : (1, 2)
    lacking = WARPSIZE*cld(size(array, cols), WARPSIZE) - size(array, cols)
    A = hcat(transpose ? Matrix(array') : array, Matrix{T}(undef, size(array, rows), lacking))
    return A
end

@inline minimum_image(s) = s - round(s)

@inline ⋅(x, y) = sum(x.*y)

function compute_tile!(forces, energies, virials, positions, L, tiles, model,
                       atoms::CUDA.CuDeviceArray{Atom,1}, ::Val{bitmask}
                       ) where {Atom, bitmask}
    @inbounds begin
        tid = CUDA.threadIdx().x
        bid = CUDA.blockIdx().x
        block_I, block_J = tiles[bid]
        is_diagonal_tile = block_J == block_I
        I = (block_I - 1)*WARPSIZE + tid
        J = (block_J - 1)*WARPSIZE + tid
        Ix3 = 3I
        Jx3 = 3J

        atom_i = atoms[I]
        atom_j = atoms[J]

        scaled_pos_i = Vec3(positions[Ix3-2], positions[Ix3-1], positions[Ix3])/L
        scaled_pos_j = Vec3(positions[Jx3-2], positions[Jx3-1], positions[Jx3])/L

        FORCES & bitmask ≠ 0 && (force_i = force_j = zeros(Vec3))
        ENERGIES & bitmask ≠ 0 && (energy_i = energy_j = 0.0f0)
        VIRIALS & bitmask ≠ 0 && (virial_i = virial_j = 0.0f0)

        for m = 1:WARPSIZE-Int(is_diagonal_tile)
            j = (tid + m - 1)%WARPSIZE + 1
            k = (tid + WARPSIZE - m - 1)%WARPSIZE + 1
            rᵥ = L*minimum_image.(scaled_pos_i - shfl_sync(CUDA.FULL_MASK, scaled_pos_j, j))
            r² = rᵥ⋅rᵥ
            E, minus_rE′ = interaction(r², model, atom_i, shfl_sync(CUDA.FULL_MASK, atom_j, j))
            if FORCES & bitmask ≠ 0
                force_ij = minus_rE′/r²*rᵥ
                force_i += force_ij
                force_j -= shfl_sync(CUDA.FULL_MASK, force_ij, k)
            end
            if ENERGIES & bitmask ≠ 0
                energy_i += E
                energy_j += CUDA.shfl_sync(CUDA.FULL_MASK, E, k)
            end
            if VIRIALS & bitmask ≠ 0
                virial_i += minus_rE′
                virial_j += CUDA.shfl_sync(CUDA.FULL_MASK, minus_rE′, k)
            end
        end

        if FORCES & bitmask ≠ 0
            CUDA.atomic_add!(pointer(forces, Ix3-2), force_i.x)
            CUDA.atomic_add!(pointer(forces, Ix3-1), force_i.y)
            CUDA.atomic_add!(pointer(forces, Ix3), force_i.z)
        end
        ENERGIES & bitmask ≠ 0 && CUDA.atomic_add!(pointer(energies, I), 0.5f0*energy_i)
        VIRIALS & bitmask ≠ 0 && CUDA.atomic_add!(pointer(virials, I), 0.5f0*virial_i)

        if !is_diagonal_tile
            if FORCES & bitmask ≠ 0
                CUDA.atomic_add!(pointer(forces, Jx3-2), force_j.x)
                CUDA.atomic_add!(pointer(forces, Jx3-1), force_j.y)
                CUDA.atomic_add!(pointer(forces, Jx3), force_j.z)
            end
            ENERGIES & bitmask ≠ 0 && CUDA.atomic_add!(pointer(energies, J), 0.5f0*energy_j)
            VIRIALS & bitmask ≠ 0 && CUDA.atomic_add!(pointer(virials, J), 0.5f0*virial_j)
        end
    end
    return nothing
end

function compute_nonbonded!(forces, energies, virials, positions, L,
                            tiles, model, atoms::CUDA.CuArray{Atom,1},
                            ::Val{bitmask}) where {Atom, bitmask}
    FORCES & bitmask ≠ 0 && (forces .= 0.0f0)
    ENERGIES & bitmask ≠ 0 && (energies .= 0.0f0)
    VIRIALS & bitmask ≠ 0 && (virials .= 0.0f0)
    CUDA.@cuda(
        threads=WARPSIZE,
        blocks=length(tiles),
        compute_tile!(forces, energies, virials, positions, L, tiles, model, atoms, Val(bitmask))
    )
end

function naively_compute_nonbonded!(forces, energies, virials, positions, L, model, atoms)
    atoms_h = Array(atoms)
    scaled_positions =map(x->x/L, Array(positions))
    N = size(positions, 2)
    forces_h = zeros(Float32, 3, N)
    energies_h = zeros(Float32, N)
    virials_h = zeros(Float32, N)
    for i = 1:N-1
        atom_i = atoms_h[i]
        scaled_position_i = scaled_positions[:,i]
        energy_i = 0.0
        virial_i = 0.0
        force_i = zeros(3)
        for j = i+1:N
            rᵥ = L*minimum_image.(scaled_position_i - scaled_positions[:,j])
            r² = rᵥ⋅rᵥ
            E, minus_E′r = interaction(r², model, atom_i, atoms_h[j])
            f_ij = minus_E′r/r²*rᵥ
            force_i += f_ij
            forces_h[:,j] -= f_ij
            energy_i += E/2
            energies_h[j] += E/2
            virial_i += minus_E′r/2
            virials_h[j] += minus_E′r/2
        end
        forces_h[:,i] += force_i
        energies_h[i] += energy_i
        virials_h[i] += virial_i
    end
    unsafe_copyto!(pointer(forces), pointer(forces_h), 3*N)
    unsafe_copyto!(pointer(energies), pointer(energies_h), N)
    unsafe_copyto!(pointer(virials), pointer(virials_h), N)
    return nothing
end

ext: jl


text:
#========================================================================================#
"""
	IN500SimpleWorld

Module IN500SimpleWorld: A Very Simple example of how to set up agents in a world.

Author: Niall Palfreyman, 26/04/22
"""
module SimpleWorld

export demo									# Externally available names of SimpleWorld
using Agents								# Required packages

#-----------------------------------------------------------------------------------------
# Module definitions:

"""
	Turtle

To construct a new agent implementation, derive it from AbstractAgent and add all necessary
attributes. The @gent macro does this automatically, also inserting the fields id, pos and vel.
speed is our own field, defining distance of travel in one step. Our use of the @agent macro
is equivalent to the following code:

	mutable struct Turtle <: AbstractAgent
		id::Int                             # Boid identity
		pos::NTuple{2,Float64}              # Position of boid
		vel::NTuple{2,Float64}              # Velocity direction of boid
		speed::Float64                      # Speed of boid
	end
"""
@agent Turtle ContinuousAgent{2} begin
	speed::Float64
end

#-----------------------------------------------------------------------------------------
# Module methods:

#-----------------------------------------------------------------------------------------
"""
	simpleworld(n)

Create a simple world containing n agents at random locations.
"""
function simpleworld( n::Int)
	# The world is continuous and 2-D, divided into 100×100 cells of width 1.0:
	world = ABM( Turtle, ContinuousSpace( (100,100), 1.0))

	# Add the n agents:
	for _ in 1:n
		add_agent!( world, Tuple(rand(1:5,2)), rand(1:9))
	end

	world
end

#-----------------------------------------------------------------------------------------
"""
	demo()

Create and run a simple world.
"""
function demo()
	# The world is continuous and 2-D, divided into 100×100 cells of width 1.0:
	world = simpleworld(5)

	println( "Number of agents before culling: ",length(allagents(world)))
	for ag in allagents(world)
		println("Agent $(ag.id) at location $(ag.pos) has speed $(ag.speed).")
	end

	for id in collect(nearby_ids( (50,50), world, 30))
		# Remove all agents within 30 units of location (50,50):
		kill_agent!( id, world)
	end

	println( "\nNumber of agents after culling: ",length(allagents(world)))
	for ag in allagents(world)
		println("Agent $(ag.id) at location $(ag.pos) has speed $(ag.speed).")
	end

end

end		# ... of module SimpleWorld


text:
using DelaySSAToolkit
using Test, SafeTestsets

@time begin
    @time @safetestset "Algorithm accuracy" begin
        include("bursty_model.jl")
    end
    @time @safetestset "Index test" begin
        include("check_index_reactions.jl")
    end
end


text:
# tools used to animate objects in the visualized scene

function (as::ArcPointsRangeSolve)(res::Vector{Float64}, x::Vector{Float64})
  res[1] = norm(x-as.x1)^2 - as.r^2
  res[2] = norm(x-as.x2)^2 - as.r^2
  if length(res) == 3
    res[3] = norm(x-as.x3)^2 - as.r^2
  end
  nothing
end

function animatearc(vc,
                    drmodel::DrawObject,
                    as::ArcPointsRangeSolve;
                    N::Int=100,
                    delaytime::Float64=0.05,
                    initrot::Rotation=Rotations.Quat(1.0,0,0,0),
                    from::Number=0,
                    to::Number=1  )
  #
  for t in linspace(from,to,N)
    am = parameterizeArcAffineMap(t, as, initrot=initrot )
    drmodel(vc, am )
    sleep(delaytime)
  end
  nothing
end

function findaxiscenter!(as::ArcPointsRangeSolve)
  d = length(as.center)
  x0 = 0.5*(as.x1+as.x2)
  r = nlsolve(as, x0)
  as.center = r.zero
  vA, vB, vC = as.x1-as.center, as.x2-as.center, as.x3-as.center
  l1, l2 = norm(as.x1-as.x2), norm(as.x2-as.x3)
  halfl0 = 0.5*norm(as.x1-as.x3)
  axt = l1 < l2 ? LinearAlgebra.cross(vA,vB) : LinearAlgebra.cross(vB,vC)
  as.axis[1:3] = axt / norm(axt)
  ta = LinearAlgebra.cross(vA,vC)
  ta ./= norm(ta)
  alph = acos(halfl0/as.r)
  if norm(ta-as.axis) < 1e-4
    #accute
    as.angle = pi - 2*alph
  else
    # oblique
    as.angle = pi + 2*alph
  end
  r.f_converged
end


text:
module TrialsTests

using BenchmarkExt
using ReTest
using StableRNGs

@testset "Trial" begin
    trial1 = BenchmarkExt.Trial(BenchmarkExt.Parameters(evals = 2))
    push!(trial1, 2, 1, 4, 5)
    push!(trial1, 21, 0, 41, 51)

    trial2 = BenchmarkExt.Trial(BenchmarkExt.Parameters(time_tolerance = 0.15))
    push!(trial2, 21, 0, 41, 51)
    push!(trial2, 2, 1, 4, 5)

    push!(trial2, 21, 0, 41, 51)
    @test length(trial2) == 3
    deleteat!(trial2, 3)
    @test length(trial1) == length(trial2) == 2
    sort!(trial2)

    @test trial1.params == BenchmarkExt.Parameters(evals = trial1.params.evals)
    @test trial2.params == BenchmarkExt.Parameters(time_tolerance = trial2.params.time_tolerance)
    @test trial1.times == trial2.times == [2.0, 21.0]
    @test trial1.gctimes == trial2.gctimes == [1.0, 0.0]
    @test trial1.memory == trial2.memory ==  [4, 41]
    @test trial1.allocs == trial2.allocs == [5, 51]

    trial2.params = trial1.params

    @test trial1 == trial2

    @test trial1[2] == push!(BenchmarkExt.Trial(BenchmarkExt.Parameters(evals = 2)), 21, 0, 41, 51)
    @test trial1[1:end] == trial1

    @test time(trial1) == time(trial2) == 2.0
    @test gctime(trial1) == gctime(trial2) == 1.0
    @test memory(trial1) == memory(trial2) == 4.0
    @test allocs(trial1) == allocs(trial2) == 5.0
    @test params(trial1) == params(trial2) == trial1.params

    # outlier trimming
    trial3 = BenchmarkExt.Trial(BenchmarkExt.Parameters(),
                                [1, 2, 3, 10, 11],
                                [1, 1, 1, 1, 1],
                                [1, 1, 1, 1, 1],
                                [1, 1, 1, 1, 1])

    trimtrial3 = rmskew(trial3)
    rmskew!(trial3)

    @test mean(trimtrial3) <= median(trimtrial3)
    @test trimtrial3 == trial3
end

@testset "TrialEstimate" begin
    rng = StableRNG(22022022)
    randtrial = BenchmarkExt.Trial(BenchmarkExt.Parameters())

    for _ in 1:40
        push!(randtrial, rand(rng, 1:20), 1, 1, 1)
    end

    while mean(randtrial) <= median(randtrial)
        push!(randtrial, rand(rng, 10:20), 1, 1, 1)
    end

    rmskew!(randtrial)

    tmin = minimum(randtrial)
    tmed = median(randtrial)
    tmean = mean(randtrial)
    tmax = maximum(randtrial)

    @test time(tmin) == time(randtrial)
    @test gctime(tmin) == gctime(randtrial)
    @test memory(tmin) == memory(tmed) == memory(tmean) == memory(tmax) == memory(randtrial)
    @test allocs(tmin) == allocs(tmed) == allocs(tmean) == allocs(tmax) == allocs(randtrial)
    @test params(tmin) == params(tmed) == params(tmean) == params(tmax) == params(randtrial)

    @test tmin <= tmed
    @test tmean <= tmed # this should be true since we called rmoutliers!(randtrial) earlier
    @test tmed <= tmax
end

@testset "TrialRatio" begin
    randrange = 1.0:0.01:10.0
    x, y = rand(randrange), rand(randrange)

    @test (ratio(x, y) == x/y) && (ratio(y, x) == y/x)
    @test (ratio(x, x) == 1.0) && (ratio(y, y) == 1.0)
    @test ratio(0.0, 0.0) == 1.0

    ta = BenchmarkExt.TrialEstimate(BenchmarkExt.Parameters(), rand(), rand(), rand(Int), rand(Int))
    tb = BenchmarkExt.TrialEstimate(BenchmarkExt.Parameters(), rand(), rand(), rand(Int), rand(Int))
    tr = ratio(ta, tb)

    @test time(tr) == ratio(time(ta), time(tb))
    @test gctime(tr) == ratio(gctime(ta), gctime(tb))
    @test memory(tr) == ratio(memory(ta), memory(tb))
    @test allocs(tr) == ratio(allocs(ta), allocs(tb))
    @test params(tr) == params(ta) == params(tb)

    @test BenchmarkExt.gcratio(ta) == ratio(gctime(ta), time(ta))
    @test BenchmarkExt.gcratio(tb) == ratio(gctime(tb), time(tb))
end

@testset "TrialJudgement" begin
    ta = BenchmarkExt.TrialEstimate(BenchmarkExt.Parameters(time_tolerance = 0.50, memory_tolerance = 0.50), 0.49, 0.0, 2, 1)
    tb = BenchmarkExt.TrialEstimate(BenchmarkExt.Parameters(time_tolerance = 0.05, memory_tolerance = 0.05), 1.00, 0.0, 1, 1)
    tr = ratio(ta, tb)
    tj_ab = judge(ta, tb)
    tj_r = judge(tr)

    @test ratio(tj_ab) == ratio(tj_r) == tr
    @test time(tj_ab) == time(tj_r) == :improvement
    @test memory(tj_ab) == memory(tj_r) == :regression
    @test tj_ab == tj_r

    tj_ab_2 = judge(ta, tb; time_tolerance = 2.0, memory_tolerance = 2.0)
    tj_r_2 = judge(tr; time_tolerance = 2.0, memory_tolerance = 2.0)

    @test tj_ab_2 == tj_r_2
    @test ratio(tj_ab_2) == ratio(tj_r_2)
    @test time(tj_ab_2) == time(tj_r_2) == :invariant
    @test memory(tj_ab_2) == memory(tj_r_2) == :invariant

    @test !(isinvariant(tj_ab))
    @test !(isinvariant(tj_r))
    @test isinvariant(tj_ab_2)
    @test isinvariant(tj_r_2)

    @test !(isinvariant(time, tj_ab))
    @test !(isinvariant(time, tj_r))
    @test isinvariant(time, tj_ab_2)
    @test isinvariant(time, tj_r_2)

    @test !(isinvariant(memory, tj_ab))
    @test !(isinvariant(memory, tj_r))
    @test isinvariant(memory, tj_ab_2)
    @test isinvariant(memory, tj_r_2)

    @test isregression(tj_ab)
    @test isregression(tj_r)
    @test !(isregression(tj_ab_2))
    @test !(isregression(tj_r_2))

    @test !(isregression(time, tj_ab))
    @test !(isregression(time, tj_r))
    @test !(isregression(time, tj_ab_2))
    @test !(isregression(time, tj_r_2))

    @test isregression(memory, tj_ab)
    @test isregression(memory, tj_r)
    @test !(isregression(memory, tj_ab_2))
    @test !(isregression(memory, tj_r_2))

    @test isimprovement(tj_ab)
    @test isimprovement(tj_r)
    @test !(isimprovement(tj_ab_2))
    @test !(isimprovement(tj_r_2))

    @test isimprovement(time, tj_ab)
    @test isimprovement(time, tj_r)
    @test !(isimprovement(time, tj_ab_2))
    @test !(isimprovement(time, tj_r_2))

    @test !(isimprovement(memory, tj_ab))
    @test !(isimprovement(memory, tj_r))
    @test !(isimprovement(memory, tj_ab_2))
    @test !(isimprovement(memory, tj_r_2))
end

@testset "Pretty printing" begin
    ta = BenchmarkExt.TrialEstimate(BenchmarkExt.Parameters(time_tolerance = 0.50, memory_tolerance = 0.50), 0.49, 0.0, 2, 1)
    tb = BenchmarkExt.TrialEstimate(BenchmarkExt.Parameters(time_tolerance = 0.05, memory_tolerance = 0.05), 1.00, 0.0, 1, 1)
    data = read(joinpath(@__DIR__, "data", "test02_pretty.txt"), String)
    pp = strip.(split(data, "\n\n\n"))

    @test BenchmarkExt.prettypercent(.3120123) == "31.20%"

    @test BenchmarkExt.prettydiff(0.0) == "-100.00%"
    @test BenchmarkExt.prettydiff(1.0) == "+0.00%"
    @test BenchmarkExt.prettydiff(2.0) == "+100.00%"

    @test BenchmarkExt.prettytime(999) == "999.000 ns"
    @test BenchmarkExt.prettytime(1000) == "1.000 μs"
    @test BenchmarkExt.prettytime(999_999) == "999.999 μs"
    @test BenchmarkExt.prettytime(1_000_000) == "1.000 ms"
    @test BenchmarkExt.prettytime(999_999_999) == "1000.000 ms"
    @test BenchmarkExt.prettytime(1_000_000_000) == "1.000 s"

    @test BenchmarkExt.prettymemory(1023) == "1023 bytes"
    @test BenchmarkExt.prettymemory(1024) == "1.00 KiB"
    @test BenchmarkExt.prettymemory(1048575) == "1024.00 KiB"
    @test BenchmarkExt.prettymemory(1048576) == "1.00 MiB"
    @test BenchmarkExt.prettymemory(1073741823) == "1024.00 MiB"
    @test BenchmarkExt.prettymemory(1073741824) == "1.00 GiB"

    @test sprint(show, "text/plain", ta) == sprint(show, ta; context=:compact => false) == pp[1]

    @test sprint(show, ta) == "TrialEstimate(0.490 ns)"
    @test sprint(
        show, ta;
        context = IOContext(
            devnull, :compact => true, :typeinfo => BenchmarkExt.TrialEstimate)
    ) == "0.490 ns"

    @test sprint(show, [ta, tb]) == "BenchmarkExt.TrialEstimate[0.490 ns, 1.000 ns]"

    trial1sample = BenchmarkExt.Trial(BenchmarkExt.Parameters(), [1], [1], [1], [1])
    @test try display(trial1sample); true catch e false end

    @static if VERSION < v"1.6-"
        @test sprint(show, "text/plain", [ta, tb]) == pp[2]
    else
        @test sprint(show, "text/plain", [ta, tb]) == pp[3]
    end

    trial = BenchmarkExt.Trial(BenchmarkExt.Parameters(), [1.0, 1.01], [0.0, 0.0], [0, 0], [0, 0])

    BenchmarkExt.set_preferences!(benchmark_output = "fancy", benchmark_histogram = "fancy")
    @test sprint(show, "text/plain", trial) == pp[4]

    BenchmarkExt.set_preferences!(benchmark_output = "classical", benchmark_histogram = "classical")
    @test sprint(show, "text/plain", trial) == pp[5]
end

end # module


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

# Devide the ErrorClassification Client, Concept, Query and Internal is possible
# with the following abstrat types.
abstract type AbstractGeneralError end
abstract type AbstractClientError <: AbstractGeneralError end
abstract type AbstractConceptError <: AbstractGeneralError end
abstract type AbstractQueryError <: AbstractGeneralError end
abstract type AbstractInternalError <: AbstractGeneralError end
abstract type AbstractServerError <: AbstractGeneralError end

# The Symbol left side is assigned to the struct which will be build bellow
# The right side of the Pair is assigned to the abstract type representing the section of error
const ERROR_STRUCTS = (
    :CLIENT_CLIENT_CLOSED=>:AbstractClientError,
    :CLIENT_SESSION_CLOSED=>:AbstractClientError,
    :CLIENT_TRANSACTION_CLOSED=>:AbstractClientError,
    :CLIENT_UNABLE_TO_CONNECT=>:AbstractClientError,
    :CLIENT_NEGATIVE_VALUE_NOT_ALLOWED=>:AbstractClientError,
    :CLIENT_MISSING_DB_NAME=>:AbstractClientError,
    :CLIENT_DB_DOES_NOT_EXIST=>:AbstractClientError,
    :CLIENT_MISSING_RESPONSE=>:AbstractClientError,
    :CLIENT_UNKNOWN_REQUEST_ID=>:AbstractClientError,
    :CLIENT_CLUSTER_NO_PRIMARY_REPLICA_YET=>:AbstractClientError,
    :CLIENT_CLUSTER_UNABLE_TO_CONNECT=>:AbstractClientError,
    :CLIENT_CLUSTER_REPLICA_NOT_PRIMARY=>:AbstractClientError,
    :CLIENT_CLUSTER_ALL_NODES_FAILED=>:AbstractClientError,
    :CONCEPT_INVALID_CONCEPT_CASTING=>:AbstractConceptError,
    :CONCEPT_MISSING_TRANSACTION=>:AbstractConceptError,
    :CONCEPT_MISSING_IID=>:AbstractConceptError,
    :CONCEPT_MISSING_LABEL=>:AbstractConceptError,
    :CONCEPT_BAD_ENCODING=>:AbstractConceptError,
    :CONCEPT_BAD_VALUE_TYPE=>:AbstractConceptError,
    :CONCEPT_BAD_ATTRIBUTE_VALUE=>:AbstractConceptError,
    :QUERY_VARIABLE_DOES_NOT_EXIST=>:AbstractQueryError,
    :QUERY_NO_EXPLANATION=>:AbstractQueryError,
    :QUERY_BAD_ANSWER_TYPE=>:AbstractQueryError,
    :QUERY_MISSING_ANSWER=>:AbstractQueryError,
    :INTERNAL_UNEXPECTED_INTERRUPTION=>:AbstractInternalError,
    :INTERNAL_ILLEGAL_STATE=>:AbstractInternalError,
    :INTERNAL_ILLEGAL_ARGUMENT=>:AbstractInternalError,
    :INTERNAL_ILLEGAL_CAST=>:AbstractInternalError,
    :GENERAL_UNKOWN_ERROR=>:AbstractInternalError,
    :GRPC_SERVER_ERROR=>:AbstractServerError)

# building the error structs according the Pairs above with some meta programming
for (T,A) in ERROR_STRUCTS
    @eval mutable struct $T <: $A
        code_prefix::String
        code_number::Int
        message_prefix::String
        message_body::String
    end
end

function Base.show(io::IO, err_struct::AbstractGeneralError)
    str = "[$(err_struct.code_prefix)$(lpad(err_struct.code_number,3,"0"))] $(err_struct.message_body)"
    print(io,str)
    return nothing
end

# All structs represent a unique number in a section of errors and the description. Both of them are
# assigned on the right side pair
const error_messages = Dict([
    #Client Error Section
    CLIENT_CLIENT_CLOSED=>(1, "The client has been closed and no further operation is allowed."),
    CLIENT_SESSION_CLOSED=>( 2, "The session has been closed and no further operation is allowed."),
    CLIENT_TRANSACTION_CLOSED=>( 3, "The transaction has been closed and no further operation is allowed."),
    CLIENT_UNABLE_TO_CONNECT=>(4, "Unable to connect to TypeDB server."),
    CLIENT_NEGATIVE_VALUE_NOT_ALLOWED=>(5, "Value cannot be less than 1, was: _error_item."),
    CLIENT_MISSING_DB_NAME=>(6,  "Database name cannot be null."),
    CLIENT_DB_DOES_NOT_EXIST=>(7,  "The database _error_item does not exist."),
    CLIENT_MISSING_RESPONSE=>(8,  "Unexpected empty response for request ID _error_item."),
    CLIENT_UNKNOWN_REQUEST_ID=>(9,  "Received a response with unknown request id _error_item."),
    CLIENT_CLUSTER_NO_PRIMARY_REPLICA_YET=>(10,  "No replica has been marked as the primary replica for latest known term _error_item."),
    CLIENT_CLUSTER_UNABLE_TO_CONNECT=>(11,  "Unable to connect to TypeDB Cluster. Attempted connecting to the cluster members, but none are available: _error_item."),
    CLIENT_CLUSTER_REPLICA_NOT_PRIMARY=>(12,  "The replica is not the primary replica."),
    CLIENT_CLUSTER_ALL_NODES_FAILED=>(13,  "Attempted connecting to all cluster members, but the following errors occurred: \n%s"),
    ### Conept Error Section
    CONCEPT_INVALID_CONCEPT_CASTING=>(1, "Invalid concept conversion from _error_item to _error_item."),
    CONCEPT_MISSING_TRANSACTION=>(2,  "Transaction cannot be null."),
    CONCEPT_MISSING_IID=>(3,  "IID cannot be null or empty."),
    CONCEPT_MISSING_LABEL=>(4,  "Label cannot be null or empty."),
    CONCEPT_BAD_ENCODING=>(5,  "The encoding _error_item was not recognised."),
    CONCEPT_BAD_VALUE_TYPE=>(6,  "The value type _error_item was not recognised."),
    CONCEPT_BAD_ATTRIBUTE_VALUE=>(7,  "The attribute value _error_item was not recognised."),
    # Query Error Section
    QUERY_VARIABLE_DOES_NOT_EXIST=>(1,  "The variable _error_item does not exist."),
    QUERY_NO_EXPLANATION=>(2,  "No explanation was found."),
    QUERY_BAD_ANSWER_TYPE=>(3,  "The answer type _error_item was not recognised."),
    QUERY_MISSING_ANSWER=>(4,  "The required field 'answer' of type _error_item was not set."),
    #Internal Error Section
    INTERNAL_UNEXPECTED_INTERRUPTION =>(1, "Unexpected thread interruption!"),
    INTERNAL_ILLEGAL_STATE =>(2, "Illegal state has been reached!"),
    INTERNAL_ILLEGAL_ARGUMENT =>(3, "Illegal argument provided: _error_item"),
    INTERNAL_ILLEGAL_CAST =>(4, "Illegal casting operation to _error_item."),
    #General Error Section included in the internal section
    GENERAL_UNKOWN_ERROR => (100, "_error_item"),
    GRPC_SERVER_ERROR => (150, "_error_item")
])

# Build the error massage according the section of errors represented by the
# abstract type of each section of errors

function call_class(typ, str, class)
    items = error_messages[class]
    return class(typ, items[1], str, items[2])
end

_build_error_messages(class::Type{<:AbstractClientError}) = call_class("CLI", "Client Error", class)
_build_error_messages(class::Type{<:AbstractConceptError}) = call_class("CON", "Concept Error", class)
_build_error_messages(class::Type{<:AbstractQueryError}) = call_class("QRY", "Query Error", class)
_build_error_messages(class::Type{<:AbstractInternalError}) = call_class("INT", "Internal Error", class)
_build_error_messages(class::Type{<:AbstractServerError}) = call_class("SRV", "Server Error", class)


text:
precision(proposal, goldstandard) =
    num_common_cons(proposal, goldstandard) / num_cons(proposal)
recall(proposal, goldstandard) =
    num_common_cons(proposal, goldstandard) / num_cons(goldstandard)
f1_measure(proposal, goldstandard) =
    2 * precision(proposal, goldstandard) * recall(proposal, goldstandard) /
       (precision(proposal, goldstandard) + recall(proposal, goldstandard))

"number of constituents of tree"
function num_cons(tree::Tree)
    M = boolean_dependency_matrix(tree)
    count(x->x, M)
end

"number of common constituents"
function num_common_cons(tree1::Tree, tree2::Tree)
    A = boolean_dependency_matrix(tree1)
    B = boolean_dependency_matrix(tree2)
    sum(a && b for (a,b) in zip(A,B))
end

"number of crossed bracket pairs"
function num_crossed_brackets(tree1, tree2)
    num_cons(tree1) - num_common_cons(tree1, tree2)
end

# t1 = tree("[0maj[0maj[0maj][1maj]][0maj[7maj][0maj]]]")
# t2 = tree("[0maj[1maj[2min][1maj]][0maj[7maj][0maj]]]")
# t3 = tree("[0maj[0maj][0maj[0maj][0maj[7maj][0maj]]]]")
#
# num_crossed_brackets(t1, t3)


text:
using
    ICNF,
    AbstractDifferentiation,
    CUDA,
    ComputationalResources,
    DataFrames,
    Distributions,
    FiniteDiff,
    FiniteDifferences,
    Flux,
    ForwardDiff,
    Optimization,
    ReverseDiff,
    MLJBase,
    SciMLBase,
    Test,
    Tracker,
    Zygote

CUDA.allowscalar() do
    include("core.jl")

    @testset "Overall" begin
        include("smoke_tests.jl")
    end
end


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

mutable struct IoK8sApiFlowcontrolV1alpha1Subject <: SwaggerModel
    group::Any # spec type: Union{ Nothing, IoK8sApiFlowcontrolV1alpha1GroupSubject } # spec name: group
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    serviceAccount::Any # spec type: Union{ Nothing, IoK8sApiFlowcontrolV1alpha1ServiceAccountSubject } # spec name: serviceAccount
    user::Any # spec type: Union{ Nothing, IoK8sApiFlowcontrolV1alpha1UserSubject } # spec name: user

    function IoK8sApiFlowcontrolV1alpha1Subject(;group=nothing, kind=nothing, serviceAccount=nothing, user=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1Subject, Symbol("group"), group)
        setfield!(o, Symbol("group"), group)
        validate_property(IoK8sApiFlowcontrolV1alpha1Subject, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiFlowcontrolV1alpha1Subject, Symbol("serviceAccount"), serviceAccount)
        setfield!(o, Symbol("serviceAccount"), serviceAccount)
        validate_property(IoK8sApiFlowcontrolV1alpha1Subject, Symbol("user"), user)
        setfield!(o, Symbol("user"), user)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1Subject

const _property_map_IoK8sApiFlowcontrolV1alpha1Subject = Dict{Symbol,Symbol}(Symbol("group")=>Symbol("group"), Symbol("kind")=>Symbol("kind"), Symbol("serviceAccount")=>Symbol("serviceAccount"), Symbol("user")=>Symbol("user"))
const _property_types_IoK8sApiFlowcontrolV1alpha1Subject = Dict{Symbol,String}(Symbol("group")=>"IoK8sApiFlowcontrolV1alpha1GroupSubject", Symbol("kind")=>"String", Symbol("serviceAccount")=>"IoK8sApiFlowcontrolV1alpha1ServiceAccountSubject", Symbol("user")=>"IoK8sApiFlowcontrolV1alpha1UserSubject")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1Subject }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1Subject))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1Subject }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1Subject[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1Subject }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1Subject[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1Subject)
    (getproperty(o, Symbol("kind")) === nothing) && (return false)
    true
end

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


text:
state = MPCCAtX(zeros(10), zeros(0), cGx = [0.0], cHx = [1.0])

@test state.x == zeros(10)
@test state.fx == nothing
@test state.gx == nothing
@test state.Hx == nothing
@test state.mu == nothing
@test state.cx == nothing
@test state.Jx == nothing

@test state.lambda == zeros(0)
@test state.current_time == nothing
@test (!(state.evals == nothing))

# On vérifie que la fonction update! fonctionne
update!(state, x = ones(10), fx = 1.0, gx = ones(10))
update!(state, lambda = ones(10), current_time = 1.0)
update!(state, Hx = ones(10,10), mu = ones(10), cx = ones(10), Jx = ones(10,10))

@test (false in (state.x .== 1.0)) == false #assez bizarre comme test...
@test state.fx == 1.0
@test (false in (state.gx .== 1.0)) == false
@test (false in (state.Hx .== 1.0)) == false
@test state.mu == ones(10)
@test state.cx == ones(10)
@test (false in (state.Jx .== 1.0)) == false
@test (false in (state.lambda .== 1.0)) == false
@test state.current_time == 1.0

reinit!(state)
@test state.x == ones(10)
@test state.fx == nothing
reinit!(state, x = zeros(10))
@test state.x == zeros(10)
@test state.fx == nothing


text:
function vech(x)
    a = triu(x)
    a = a[find(a)]
end


text:
"""
    Pkn_NGG_arb(k, n, β, σ)

Compute the prior probability mass to obtain k clusters through direct computation with arbitrary precision. Not the most stable version.

# Examples
```julia-repl
julia> GibbsTypePriors.Pkn_NGG_arb(50, 100, 0.5, 0.2)
[2.14762363546438648278532704280127812426643249551735595554664760044538213621775301015485328374549413897540160511465187297274194698946624496309582378995095872169826130349407122713909436452654024908060179841056741955863382744781148806296586030164131604285251938875060389310445972984208190921903755337680167605513750214014651799055934001372969787136894687051060898423210469539540879653779017239230826328548585098581761936756315153017814207498994142860020958294285847670933928577444047966837019301853075052805215643829733023817130541787687113481647039349350666819737212357118184059914081469496628351252560981171967614071257194538987858556850738404955654861591001195199137908090799330617531307554849423121028628550128016847465035159250895414945842092857092939501917174224657462412244406519738404915932075136518078258780668230018325226513272503217141555634764038959108448484312201406445495467152388854367276670509646840248971542665576196296601294567533128012684921641261767619308293743726390603225250672867733057200093764412465619355633068601123653038002743467832229736915692389253166213741137134870700536044628527732021542164522846984454807319649387483162094391527262191665810862006754021341268115958393216764683235046467243249773739254152551430288889983085214706056491738375963632435912736357674313108409758462287962329443736179207688980017178194117395394271454788659703559151681296053952789173258126586683260147689572558384333967140709303631030995994161910062213092853666415615122632e-19 +/- 4.26e-1498]
```
"""
function Pkn_NGG_arb(k, n, β, σ)
    if k>n || k==0
        return arb_0
    else
        σ_arb = RR(σ)
        Vnk_NGG(n, k, β, σ) // σ_arb^k * Cnk(n, k, σ)
    end
end

"""
    Pkn_NGG_robust_in(k, n, β, σ)

Compute the prior probability mass to obtain k clusters, either through direct computation or using the recursive computation for the Cnk, with arbitrary precision. Stable for large values of k.

# Examples
```julia-repl
julia> GibbsTypePriors.Pkn_NGG_robust_in(50, 100, 0.5, 0.2)
[2.14762363546438648278532704280127812426643249551735595554664760044538213621775301015485328374549413897540160511465187297274194698946624496309582378995095872169826130349407122713909436452654024908060179841056741955863382744781148806296586030164131604285251938875060389310445972984208190921903755337680167605513750214014651799055934001372969787136894687051060898423210469539540879653779017239230826328548585098581761936756315153017814207498994142860020958294285847670933928577444047966837019301853075052805215643829733023817130541787687113481647039349350666819737212357118184059914081469496628351252560981171967614071257194538987858556850738404955654861591001195199137908090799330617531307554849423121028628550128016847465035159250895414945842092857092939501917174224657462412244406519738404915932075136518078258780668230018325226513272503217141555634764038959108448484312201406445495467152388854367276670509646840248971542665576196296601294567533128012684921641261767619308293743726390603225250672867733057200093764412465619355633068601123653038002743467832229736915692389253166213741137134870700536044628527732021542164522846984454807319649387483162094391527262191665810862006754021341268115958393216764683235046467243249773739254152551430288889983085214706056491738375963632435912736357674313108409758462287962329443736179207688980017178194117395394271454788659703559151681296053952789173258126586683260147689572558384333967140709303631030995994161910062213092853666415615122632e-19 +/- 4.26e-1498]
```
"""
function Pkn_NGG_robust_in(k, n, β, σ; verbose = false)
    if k>n || k==0
        return arb_0
    else
        σ_arb = RR(σ)
        Vnk_NGG_rec(n, k, β, σ) // σ_arb^k * Cnk_robust(n, k, σ; verbose = verbose)
    end
end

"""
    Pkn_NGG_robust(k, n, β, σ)

Compute the prior probability mass to obtain k clusters, either through direct computation or using the recursive computation for the Cnk, with arbitrary precision. Stable for large values of k. If this became unstable for small values of k, the function will print a warning and resort to using an approximation scheme based on the predictive distribution of Gibbs-type processes (Bystrova, 2020)

References:
- Bystrova, D, Kon Kam King, G, Deslandes, F, Arbel, J. Approximating the clusters' prior distribution in Bayesian nonparametric models, 2020 (in preparation).

# Examples
```julia-repl
julia> GibbsTypePriors.Pkn_NGG_robust(1, 100, 1.2, 0.6)
[1.9296391903064605384654157430472418076148726737727875230129641023067313090166540244753747312457129183657025493616396554177981937762780521371044832530798250427147509991595690117491772034391001735237273292732682815399768750800178863279517168277700486277214486303628778846185502064952747248507331917928973821847721861122917899722304242486554537125500815028027048749342590041556047806174326663544410878820833633145820033554401795667161239326297269828805842051466014698678958849980095944598411211901947299487314690595562677257989440632872268784586850284818761100358647158324006198078594045133536472908933706208286326023740694795701448737505642410975566985814946675218483636408975980085854487052384785613476386696535340182246327012856942430959614307701381885382101908730464363930653564743180041556516459003976577914678253451116729865709223986356119211820112193877506211938987401933281878273977621027890553237363870948374184877771847370494105304641393041184948881280139977645736656738738876663695255598855102890090929354876451974185700013265883795024142293167970426325083355312585087826163789033506198771766732204008615074408731001784408725973789012051139760486038536969512748787960733295275525464396278674073554207828665628723052100436461108870416804992391459295951069501611168607219137361112920188489260042042387173028715544129182561207314706608564475318060162678254188901156369626393351965452285428586670794557016665659781242771755478946260451052238e-6 +/- 5.76e-1451]
```
"""
function Pkn_NGG_robust(k, n, β, σ; verbose = false)
    res = Pkn_NGG_robust_in(k, n, β, σ; verbose = verbose)
    if has_reasonable_precision(res)
        return res
    else
        @warn "Direct computation not accurate enough, resorting to approximation"
        start_Pkn_val = res
        i = 1
        start_k_ind = undef
        while !has_reasonable_precision(start_Pkn_val)
            if 2^i*k>n
                error("There seem to be a numerical problem with computing Pkn even for large values of k")
            end
            start_k_ind = 2^i*k
            start_Pkn_val = Pkn_NGG_robust_in(start_k_ind, n, β, σ; verbose = false)#verbose = false bccause this can be very informative
            i += 1
        end
        if verbose
            println("k=$k, start_k_ind=$start_k_ind")
        end
        return Pkn_NGG_approx(k, n, β, σ, start_k_ind, start_Pkn_val)
    end
end

"""
    log_βnk(β, n, k, σ)

log-factor involved in the expansion of the predictive distribution (Bystrova, 2020).

References:
- Bystrova, D, Kon Kam King, G, Deslandes, F, Arbel, J. Approximating the clusters' prior distribution in Bayesian nonparametric models, 2020 (in preparation).
```
"""
function log_βnk(β, n, k, σ)
    return log(β) + log(n) - 1/σ*log(k)
end
"""
    βnk(β, n, k, σ)

    factor involved in the expansion of the predictive distribution (Bystrova, 2020).

    References:
    - Bystrova, D, Kon Kam King, G, Deslandes, F, Arbel, J. Approximating the clusters' prior distribution in Bayesian nonparametric models, 2020 (in preparation).
```
"""
βnk(β, n, k, σ) = exp(log_βnk(β, n, k, σ))

"""
    logxk(n, k, β, σ)

log of the ratio p_{n+1,k+1}/p_{n+1,k} (Bystrova, 2020).

    References:
    - Bystrova, D, Kon Kam King, G, Deslandes, F, Arbel, J. Approximating the clusters' prior distribution in Bayesian nonparametric models, 2020 (in preparation).
```
"""
function logxk(n, k, β, σ)
    if n==1
        return log(k*σ ) + log(Cnk(n, k+1, σ)) - log(σ) - log(Cnk(n, k, σ))
    else
       return log(k*σ + βnk(β, n-1, k, σ)) + log(Cnk(n, k+1, σ)) - log(σ) - log(Cnk(n, k, σ))
    end
end

"""
    logxk_py(n, k, β, σ)

log of the ratio p_{n+1,k+1}/p_{n+1,k} for preddictive approximation for Poisson-Dirichlet process (Bystrova, 2020).

    References:
    - Bystrova, D, Kon Kam King, G, Deslandes, F, Arbel, J. Approximating the clusters' prior distribution in Bayesian nonparametric models, 2020 (in preparation).
```
"""

function logxk_py(n, k, β, σ)
       return log(k*σ + β) + log(Cnk(n, k+1, σ)) - log(σ) - log(Cnk(n, k, σ))
end

"""
    Pkn_NGG_approx(k, n, β, σ, start_k_ind, start_Pkn_val)

Compute the prior probability mass to obtain k clusters using an approximation scheme based on the predictive distribution of Gibbs-type processes (Bystrova, 2020). This starts from a hopefully exact value for large k and recursively computes values for smaller k.

References:
- Bystrova, D, Kon Kam King, G, Deslandes, F, Arbel, J. Approximating the clusters' prior distribution in Bayesian nonparametric models, 2020 (in preparation).

# Examples
```julia-repl
julia> GibbsTypePriors.Pkn_PY(3, 20, 1.2, 0.6, 4, 0.5)
```
"""
function Pkn_NGG_approx(k, n, β, σ, start_k_ind, start_Pkn_val)
    if k==start_k_ind
        return start_Pkn_val
    else
        logPkn = log(start_Pkn_val)
        for ki in (start_k_ind-1):-1:k
            # println(logPkn)
            logPkn = logPkn - logxk(n, ki, β, σ)
        end
    #     return P1kn
    #     [logprod_xi(n, j, β, σ) for j in 1:(n-1)]
        return exp(logPkn)
    end
end

"""
    Pkn_NGG_approx(k, n, β, σ)

Compute the prior probability mass to obtain k clusters using an approximation scheme based on the predictive distribution of Gibbs-type processes (Bystrova, 2020). This starts from a hopefully exact value for k = n and recursively computes values for smaller k.

References:
- Bystrova, D, Kon Kam King, G, Deslandes, F, Arbel, J. Approximating the clusters' prior distribution in Bayesian nonparametric models, 2020 (in preparation).

# Examples
```julia-repl
julia> GibbsTypePriors.Pkn_NGG_approx(1, 20, 1.2, 0.6)
[2.7009399464328485036455096714993560085838162743940856587071430940695692450282156669202521580688533184574332460459428922676018622596226389639876939210136324011290876441564860327431322852335827166684153681574741014245777351580276219866801498262214575832111143647387542038414573675796133766555477716597201505157807253229752694991002241240470559749206273197429789323911087400986000916736286072849665709327518651195311810713486088342136293790842083300605780362633828289343957046209090018048659874748297475822372662173612817812499579367976643357372413282552630908029675874998289474675856473298883544580045158165555501162265151669548888200328109680536750079782964939851007608028946482880684475495853392663598479373751535177124555359076603669329448695955740396252130711942753713445032238102287762219833145535097175721479527568503456921692682357547219562637788922234866812786408688403987111183409338532478791363489211597726862086464685845551045734603347887381596652871622921982618325894422172400624261984839422035526236600222545483736019629179100904287763109619430949410201339488837712451860950411175609450393575033561080575062498453159378887705857471003486442931337325289901768958343433009545255517916504968781730504152424300613835412129901208365234370965057437662018159597238815095679061796910391633816437061299873914601373715585590808066739330963287595137055126014533391777974177238039468992331265855933979299517979422838595919693752234936218541025748729993033182387854702671683e-5 +/- 4.66e-1477]
```
"""
Pkn_NGG_approx(k, n, β, σ) = Pkn_NGG_approx(k, n, β, σ, Int64(floor(n/2)), Pkn_NGG_arb(Int64(floor(n/2)), n, β, σ))

"""
    Pkn_NGG_approx(k, n, β, σ, start_k_ind, start_Pkn_val)

Compute the prior probability mass to obtain k clusters for all k between 1 and start_Pkn_val using an approximation scheme based on the predictive distribution of Gibbs-type processes (Bystrova, 2020). This starts from a hopefully exact value for k = n and recursively computes values for smaller k.

References:
- Bystrova, D, Kon Kam King, G, Deslandes, F, Arbel, J. Approximating the clusters' prior distribution in Bayesian nonparametric models, 2020 (in preparation).

# Examples
```julia-repl
julia> GibbsTypePriors.Pkn_NGG_approx(20, 1.2, 0.6)
```
"""
function Pkn_NGG_approx(n, β, σ, start_k_ind, start_Pkn_val)
    logP1kn = Array{arb}(undef, start_k_ind)
    logP1kn[start_k_ind] = start_Pkn_val |> log
    for k in (start_k_ind-1):-1:1
        logP1kn[k] = logP1kn[k+1] - logxk(n, k, β, σ)
    end
    if start_k_ind == n
        return  exp.(logP1kn)./(sum(exp.(logP1kn)))
    else
      return exp.(logP1kn)
    end
end

"""
    Pkn_NGG_approx(n, β, σ)

Compute the prior probability mass to obtain k clusters for all k between 1 and n using an approximation scheme based on the predictive distribution of Gibbs-type processes (Bystrova, 2020). This starts from a hopefully exact value for k = n and recursively computes values for smaller k.

References:
- Bystrova, D, Kon Kam King, G, Deslandes, F, Arbel, J. Approximating the clusters' prior distribution in Bayesian nonparametric models, 2020 (in preparation).

# Examples
```julia-repl
julia> GibbsTypePriors.Pkn_NGG_approx(20, 1.2, 0.6)
```
"""
Pkn_NGG_approx(n, β, σ) = Pkn_NGG_approx(n, β, σ, n, Pkn_NGG_arb(n, n, β, σ))

function Pkn_NGG_robust(n, β, σ; verbose = false)
    P1n = Array{arb}(undef, n)
    int_n_half::Int64 = floor(n/2)
    P1n[int_n_half:n] = Pkn_NGG_arb.(int_n_half:n, n, β, σ)
    Pkn_val = P1n[int_n_half]
    if !has_reasonable_precision(Pkn_val)
        error("There seem to be a numerical problem with computing Pkn even for k=n/2")
    else
        k = int_n_half
        while k ≥ 1 && has_reasonable_precision(Pkn_val)
            P1n[k] = Pkn_val
            k -= 1
            Pkn_val = Pkn_NGG_arb(k, n, β, σ)
        end
        if verbose
            println("k=$k, n=$n")
        end
        if k > 1
            @warn "Direct computation not accurate enough, resorting to approximation"
            P1n[1:(k+1)] = Pkn_NGG_approx(n, β, σ, k+1, P1n[k+1])
        end
        return P1n
    end
end

"""
    Pkn_NGG(k, n, β, σ)

Compute the prior probability mass to obtain k clusters, either through direct computation or using the recursive computation for the Cnk, with arbitrary precision. Stable for large values of k, the function may become unstable for small values of k. If this were the case, the function would print a warning and resort to using an approximation scheme based on the predictive distribution of Gibbs-type processes (Bystrova, 2020)

References:
- Bystrova, D, Kon Kam King, G, Deslandes, F, Arbel, J. Approximating the clusters' prior distribution in Bayesian nonparametric models, 2020 (in preparation).

# Examples
```julia-repl
julia> GibbsTypePriors.Pkn_NGG(10, 100, 1.2, 0.6)
```
"""
Pkn_NGG(k, n, β, σ) = Pkn_NGG_robust(k, n, β, σ; verbose = false) |> Float64

"""
    Pkn_NGG(n, β, σ)

Compute the prior probability mass to obtain k clusters over k=1 to k=n, either through direct computation or using the recursive computation for the Cnk, with arbitrary precision. Stable for large values of k. If this became unstable for small values of k, the function will print a warning and resort to using an approximation scheme based on the predictive distribution of Gibbs-type processes (Bystrova, 2020)

References:
- Bystrova, D, Kon Kam King, G, Deslandes, F, Arbel, J. Approximating the clusters' prior distribution in Bayesian nonparametric models, 2020 (in preparation).

# Examples
```julia-repl
julia> GibbsTypePriors.Pkn_NGG(100, 1.2, 0.6)
```
"""
#Pkn_NGG(n, β, σ) = Pkn_NGG_robust(k,n, β, σ; verbose = false) |> Float64
Pkn_NGG(n, β, σ) = convert(Array{Float64,1},Pkn_NGG_robust.(1:n, n, β, σ; verbose = false))

"""
    Pkn_NGG_pred_approx(n, β, σ)
Compute the approximation for prior probability mass to obtain k clusters over k=1 to k=n, either using approximation scheme based on the predictive distribution of Gibbs-type processes (Bystrova, 2020) and through direct computation or using the recursive computation for the Cnk, with arbitrary precision.
References:
- Bystrova, D, Kon Kam King, G, Deslandes, F, Arbel, J. Approximating the clusters' prior distribution in Bayesian nonparametric models, 2020 (in preparation).
# Examples
```julia-repl
julia> GibbsTypePriors.Pkn_NGG_pred_approx(100, 1.2, 0.6)
```
"""

function Pkn_NGG_approx_full(n, β, σ, f)
        Axnk = Array{arb}(undef, n-1)
        Axnk[1:n-1] = f.(n, 1:n-1, β, σ)
        Sum_xn= exp(f(n, 1, β, σ))
        for i in (3:n)
             Sum_xn= Sum_xn + exp(sum(Axnk[1:i-1]))
        end
        P1n = Array{arb}(undef, n)
        P1n[1] = exp(- log(1 +Sum_xn))
        for k in 2:(n)
              P1n[k] = exp(log(P1n[k-1]) + Axnk[k-1])
        end
        return P1n
end

Pkn_NGG_pred_approx(n, β, σ) = convert(Array{Float64,1}, Pkn_NGG_approx_full(n,  β, σ, logxk))

Pkn_PY_pred_approx(n, β, σ) = convert(Array{Float64,1}, Pkn_NGG_approx_full(n,  β, σ, logxk_py))

function Pkn_2PD_arb(k::N, n::N, θ::T, σ::T) where {T<:Number, N<:Integer}
    σ_arb = RR(σ)
    Vnk_2PD(n, k, θ, σ) // σ_arb^k * Cnk(n, k, σ)
end

"""
    Pkn_2PD(n, β, σ)

Compute the prior probability mass to obtain k clusters for a 2-parameter Poisson-Dirichlet process (also known as Pitman-Yor).

# Examples
```julia-repl
julia> Pkn_2PD(10, 100, 1.2, 0.6)
0.00787749555803404
```
"""
Pkn_2PD(k::N, n::N, θ::T, σ::T) where {T<:Number, N<:Integer} = convert(Float64, Pkn_2PD_arb(k, n, θ, σ))
Pkn_PY_arb = Pkn_2PD_arb
"""
    Pkn_PY(n, θ, σ)

Compute the prior probability mass to obtain k clusters for a 2-parameter Poisson-Dirichlet process (also known as Pitman-Yor).

# Examples
```julia-repl
julia> Pkn_PY(100, 1.2, 0.6)
```
"""
Pkn_PY(n, θ, σ) =  convert(Array{Float64,1}, Pkn_PY_arb.(1:n, n, θ, σ))

function Pkn_Dirichlet_arb(k, n, θ)
    θ_arb = RR(θ)
    return θ_arb^k // risingfac(θ_arb,n) * unsigned_Stirling1(n,k)
end

"""
    Pkn_Dirichlet(k, n, θ)

Compute the prior probability mass to obtain k clusters for a Dirichlet process.

# Examples
```julia-repl
julia> Pkn_Dirichlet(7, 10, 1.5)
```
"""
Pkn_Dirichlet(k, n, θ) = convert(Float64, Pkn_Dirichlet_arb(k, n, θ))

function Pkn_Dirichlet_mult_arb(k, n, H,  θ)
    if (k<1)||(k>n)
        return 0.0
    else
        H_arb = RR(H)
        θ_arb = RR(θ)
        return (fac(H_arb)//(fac(H_arb - k)*risingfac(θ_arb,n))) *((-1)^k )* Cnk(n, k,-θ_arb//H_arb)
    end
end

"""
    Pkn_Dirichlet_mult(k, n, H, θ)
Compute the prior probability mass to obtain k clusters for a Dirichlet multinomial process.
# Examples
```julia-repl
julia> Pkn_Dirichlet_mult(7, 10,100, 1.5)
0.0100815973559968
```
"""

Pkn_Dirichlet_mult(k, n, H, θ) = convert(Float64, Pkn_Dirichlet_mult_arb(k, n,H,θ))

function Pkn_PDM_arb(k::N, n::N,H::N, θ::T, σ::T) where {T<:Number, N<:Integer}
    H_arb = RR(H)
    θ_arb = RR(θ)
    σ_arb= RR(σ)
    coef= fac(H_arb)// fac(H_arb -k)
    return coef* sum([ (Vnk_2PD(n, i, θ, σ) * unsigned_Stirling2(i,k)* Cnk(n, i, σ))// ((H_arb* σ_arb)^i)  for i in k:n])
  end

  """
      Pkn_PYM(k, n, H, θ)
  Compute the prior probability mass to obtain k clusters for a Pitma-Yor multinomial process.
  # Examples
  ```julia-repl
  julia> Pkn_PYM(7, 10, 100,  1.5, 0.2)
  0.05547962896578268
  ```
  """

Pkn_PYM(k, n, H, θ, σ) = Pkn_PDM_arb(k,n,H,θ, σ) |> Float64

function Pkn_NGGM_arb(k::N, n::N,H::N,  β::T, σ::T) where {T<:Number, N<:Integer}
    H_arb = RR(H)
    return (fac(H_arb) // (H_arb^k * fac(H_arb-k))) * sum([ (1//(H_arb^i)) * unsigned_Stirling2(i+k,k)* Pkn_NGG_arb(k+i, n,  β, σ)  for i in 0:(n-k)])
end

"""
    Pkn_NGGM(k, n, H, β, σ)
Compute the prior probability mass to obtain k clusters for a normalized generalized gamma multinomial process.
Computationally extensive for multiple k.
# Examples
```julia-repl
julia> Pkn_NGGM(7, 10, 100,  1.0, 0.25)
0.005029270051158085
```
"""
Pkn_NGGM(k, n, H, β, σ) = convert(Float64, Pkn_NGGM_arb(k, n, H, β, σ))

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
@testset "1000.minimum-cost-to-merge-stones.jl" begin
    @test merge_stones([3,2,4,1], 3) == -1
    @test merge_stones([3,5,1,2,6], 3) == 25
    @test merge_stones([3,2,4,1], 2) == 20
end


text:
using CompScienceMeshes
using Test

function isoriented(m::Mesh)

    @assert dimension(m) >= 0
    edges = skeleton(m, dimension(m)-1)

    D = connectivity(edges, m)
    S = (sum(D,1) .== 0)
    return all(S)
end

fn = joinpath(dirname(@__FILE__),"assets","cube2.in")
m = readmesh(fn)
@test isoriented(m)

fn = joinpath(dirname(@__FILE__),"assets","sphere2.in")
m = readmesh(fn)
@test isoriented(m)


text:
#=
# Uses the elemental order to determine the nodal order.
# Order within elements is lexicographic. Other options may be added.
=#
export reorder_grid_element_first!

function reorder_grid_element_first!(grid, porder, elorder)
    dim = size(grid.allnodes,1);
    new_node_order= element_first_order_map(grid, porder, elorder, dim, true);
    reorder_grid_nodes!(grid, new_node_order);
    return grid;
end

function element_first_order_map(grid, porder, elorder, dim, invert = true)
    nnodes = size(grid.allnodes,2);
    eforder = zeros(Int, nnodes); # The ordering
    node_done = zeros(Int, nnodes);

    # Loop over the elements in elorder, filling in the eforder as we go
    l2g = grid.loc2glb;
    efind = 0;
    for ei=elorder
        for ni=1:size(grid.loc2glb,1)
            # node's old index
            old_ind = grid.loc2glb[ni,ei];
            if node_done[old_ind] == 0 # not done yet
                # node's new index
                efind = efind+1;

                if invert
                    eforder[old_ind] = efind;
                else
                    eforder[efind] = old_ind;
                end
                node_done[old_ind] = efind;

            else # already mapped

            end

        end

        # # Add interior first
        # for k=1:porder+1
        #     if dim > 1
        #         for j = 1:porder+1
        #             if dim > 2
        #                 for i = 1:porder+1
        #                     # 3D
        #                     lind = i + (porder+1)*((j-1) + (porder+1)*(k-1)); # local index
        #                     gind = l2g[lind, ei]; # global index
        #                     if !node_done[gind]
        #                         efind = efind+1;
        #                         if invert
        #                             eforder[gind] = efind;
        #                         else
        #                             eforder[efind] = gind;
        #                         end
        #                         node_done[gind] = true;
        #                     end
        #                 end # i loop

        #             else # 2D
        #                 lind = j + (porder+1)*(k-1); # local index
        #                 gind = l2g[lind, ei]; # global index
        #                 if !node_done[gind]
        #                     efind = efind+1;
        #                     if invert
        #                         eforder[gind] = efind;
        #                     else
        #                         eforder[efind] = gind;
        #                     end
        #                     node_done[gind] = true;
        #                 end
        #             end

        #         end # j loop

        #     else # 1D
        #         lind = k; # local index
        #         gind = l2g[lind, ei]; # global index
        #         if !node_done[gind]
        #             efind = efind+1;
        #             if invert
        #                 eforder[gind] = efind;
        #             else
        #                 eforder[efind] = gind;
        #             end
        #             node_done[gind] = true;
        #         end
        #     end

        # end # k loop
    end # elements

    return eforder;
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
using COSMO, Random, Test, LinearAlgebra, SparseArrays, Random
rng = Random.MersenneTwister(144545)

# We create a test problem with 4 convex sets (2 of them decomposable), 1 zero set in the middle
@testset "Decomposition with PSDTriangle" begin

# define convex set
C = [COSMO.PsdCone(16); COSMO.ZeroSet(2); COSMO.PsdCone(16); COSMO.PsdCone(9)]

m = 1

# CONSTRAINT 1
A1 = rand(rng, 4, 4)
A1 = 0.5 * (A1 + A1')
# eliminate certain entries to make matrix chordal
A1[1, 3] = A1[1, 4] = A1[3, 1] = A1[4, 1] = 0
a1 = vec(A1)
# find solution matrix S that is pos def
S1 = A1 + (eigvals(A1)[1] + 1) * Matrix(1.0I, 4, 4)

# CONSTRAINT 2
a2 = rand(rng, 2, 1)
s2 = [0; 0]

# CONSTRAINT 3
A3 = rand(rng, 4, 4)
A3 = 0.5 * (A3 + A3')
# eliminate certain entries to make matrix chordal
A3[2, 4] = A3[4, 2] =  A3[1, 3]  = A3[3, 1] = A3[2, 3] = A3[3, 2] = 0
a3 = vec(A3)
S3 = A3 + (eigvals(A3)[1] + 1) * Matrix(1.0I, 4, 4)

# CONTRAINT 4 (dense and undecomposable)
A4 = rand(rng, 3, 3)
A4 = 0.5 * (A4 + A4')
a4 = vec(A4)
S4 = A4 + (eigvals(A4)[1] + 1) * Matrix(1.0I, 3, 3)

# Find b such that problem is feasible
x = rand(rng, 1)
A = [a1; a2; a3; a4]
s = [vec(S1); s2; vec(S3); vec(S4)]
b = A * x + s

# Find q such that problem is feasible
Y1 = generate_pos_def_matrix(rng, 4,  0.1, 1.)
y2 = rand(rng, 2)
Y3 = generate_pos_def_matrix(rng, 4,  0.1, 1.)
Y4 = generate_pos_def_matrix(rng, 3,  0.1, 1.)
y = [vec(Y1); y2; vec(Y3); vec(Y4)]

P = sparse(zeros(1, 1))
q = -P * x - A' * y

# --------------------------
# CONFIG 1: PsdSquare + no decomposition
# --------------------------

model = COSMO.Model()
settings = COSMO.Settings(decompose = false)
COSMO.set!(model, P, q, A, b, C, settings)
res1 = COSMO.optimize!(model);
# Y1_sol1 = reshape(res1.y[1:16], 4, 4)
# Y2_sol1 = reshape(res1.y[19:34], 4, 4)
S1_sol1 = reshape(res1.s[1:16], 4, 4)
S2_sol1 = reshape(res1.s[19:34], 4, 4)

# --------------------------
# CONFIG 2: PsdSquare +  decomposition
# --------------------------

model = COSMO.Model()
settings = COSMO.Settings(decompose = true, compact_transformation = false)
 COSMO.set!(model, P, q, A, b, C, settings)
res2 = COSMO.optimize!(model);
S1_sol2 = reshape(res2.s[1:16], 4, 4)
S2_sol2 = reshape(res2.s[19:34], 4, 4)

# --------------------------
# Only keep upper triangular data of the problem
# --------------------------

# Notice that we have to rescale the off-diagonal entries of s by sqrt(2) since our projection assumes
# that
Ct = [COSMO.PsdConeTriangle(10); COSMO.ZeroSet(2); COSMO.PsdConeTriangle(10); COSMO.PsdConeTriangle(6)]

a1t = zeros(10)
COSMO.extract_upper_triangle!(A1, a1t, sqrt(2))
a3t = zeros(10)
COSMO.extract_upper_triangle!(A3, a3t, sqrt(2))
a4t = zeros(6)
COSMO.extract_upper_triangle!(A4, a4t, sqrt(2))
At = [a1t; a2; a3t; a4t]

B1 = reshape(b[1:16], 4, 4)
B3 = reshape(b[19:34], 4, 4)
B4 = reshape(b[35:end], 3, 3)
b1t = zeros(10)
b3t = zeros(10)
b4t = zeros(6)
COSMO.extract_upper_triangle!(B1, b1t, sqrt(2))
COSMO.extract_upper_triangle!(B3, b3t, sqrt(2))
COSMO.extract_upper_triangle!(B4, b4t, sqrt(2))
bt = [b1t; b[17:18]; b3t; b4t]

# --------------------------
# CONFIG 3: PsdTriangle + no decomposition
# --------------------------

model = COSMO.Model()
settings = COSMO.Settings(decompose = false)
COSMO.set!(model, P, q, At, bt, Ct, settings)
res3 = COSMO.optimize!(model);
S1_sol3 = zeros(4, 4)
COSMO.populate_upper_triangle!(S1_sol3, res3.s[1:10], 1. / sqrt(2))
S1_sol3 = Symmetric(S1_sol3)
S2_sol3 = zeros(4, 4)
COSMO.populate_upper_triangle!(S2_sol3, res3.s[13:22], 1. / sqrt(2))
S2_sol3 = Symmetric(S2_sol3)

# --------------------------
# CONFIG 4: PsdTriangle +  decomposition
# --------------------------

model = COSMO.Model()
settings = COSMO.Settings(decompose = true)
COSMO.set!(model, P, q, At, bt, Ct, settings)
res4 = COSMO.optimize!(model);
S1_sol4 = zeros(4, 4)
COSMO.populate_upper_triangle!(S1_sol4, res4.s[1:10], 1. / sqrt(2))
S1_sol4 = Symmetric(S1_sol4)
S2_sol4 = zeros(4, 4)
COSMO.populate_upper_triangle!(S2_sol4, res4.s[13:22], 1. / sqrt(2))
S2_sol4 = Symmetric(S2_sol4)
Y1_sol4 = matrixify(res4.y[1:10])
Y2_sol4 = matrixify(res4.y[13:22])
@testset "Check correct solution" begin
  @test abs(res1.obj_val - res2.obj_val) < 1e-4
  @test abs(res1.obj_val - res3.obj_val) < 1e-4
  @test abs(res2.obj_val - res4.obj_val) < 1e-4
  @test norm(S1_sol1 - S1_sol2, Inf) < 1e-5
  @test norm(S1_sol2 - S1_sol3, Inf) < 1e-5
  @test norm(S1_sol3 - S1_sol4, Inf) < 1e-5
  @test norm(S2_sol1 - S2_sol2, Inf) < 1e-5
  @test norm(S2_sol2 - S2_sol3, Inf) < 1e-5
  @test norm(S2_sol3 - S2_sol4, Inf) < 1e-5

end

# Check that the correct y is reassembled after the decomposition, i.e. the correct dual error holds and that Y is positive semidefinite
pattern2 = ones(9, 9)
pattern2[1, 5:9] .= 0
pattern2[2, 6:9] .= 0
pattern2[3, 6:9] .= 0
pattern2[4, 9] = 0
pattern2 = Matrix(Symmetric(pattern2, :U))

P, q, At, bt, Ct = feasible_sdp_with_pattern(rng, pattern2)

settings = COSMO.Settings(decompose = false, verbose =true);
model = COSMO.Model()
COSMO.set!(model, P, q, At, bt, Ct, settings)
res5 = COSMO.optimize!(model);

settings = COSMO.Settings(decompose = true, complete_dual = true, verbose = true, compact_transformation = true);
model = COSMO.Model()
COSMO.set!(model, P, q, At, bt, Ct, settings)
res6 = COSMO.optimize!(model);

settings = COSMO.Settings(decompose = true, complete_dual = true, compact_transformation = false);
model = COSMO.Model()
COSMO.set!(model, P, q, At, bt, Ct, settings)
res7 = COSMO.optimize!(model);

Y_5 = matrixify(res5.y)
Y_6 = matrixify(res6.y)
Y_7 = matrixify(res7.y)
@testset "Check correct assembly of dual variable y" begin
  @test minimum(eigen(Y_5).values) > -1e-6
  #@test minimum(eigen(Y_6).values) > -1e-6 #TODO: Check why this errors
  #@test minimum(eigen(Y_7).values) > -1e-6
  @test norm(At' * res5.y + q, Inf) < 1e-3
  @test norm(At' * res6.y + q, Inf) < 1e-3
  @test norm(At' * res7.y + q, Inf) < 1e-3
end
end

nothing


text:
using Test
import IJulia, JSON

@testset "installkernel" begin
    let kspec = IJulia.installkernel("ijuliatest", "-O3", "-p2",
                    env=Dict("FOO"=>"yes"), specname="Yef1rLr4kXKxq9rbEh3m")
        try
            @test dirname(kspec) == IJulia.kerneldir()
            @test isfile(joinpath(kspec, "kernel.json"))
            @test isfile(joinpath(kspec, "logo-32x32.png"))
            @test isfile(joinpath(kspec, "logo-64x64.png"))
            let k = open(JSON.parse, joinpath(kspec, "kernel.json"))
                debugdesc = ccall(:jl_is_debugbuild,Cint,())==1 ? "-debug" : ""
                @test k["display_name"] == "ijuliatest" * " " * Base.VERSION_STRING * debugdesc
                @test k["argv"][end] == "{connection_file}"
                @test k["argv"][end-3:end-2] == ["-O3", "-p2"]
                @test k["language"] == "julia"
                @test k["env"]["FOO"] == "yes"
            end
        finally
            rm(kspec, force=true, recursive=true)
        end
    end
end


text:
# 0 Job version means that no jobs have ran yet.
"Returns the version found in the .metadata.jld2 if it exists. Otherwise 0."
main_job_version(dir::AbstractString) = job_versions(dir)[end]
main_job_version(job::Job) = main_job_version(abspath(job))

function job_versions(dir::AbstractString)
    verdir = joinpath(dir, VERSION_DIR_NAME)
    if ispath(verdir)
        versions = parse.(Int, readdir(verdir))
        push!(versions, versions[end]+1)
    else
        versions = Int[0]
    end

    return sort(unique(versions))
end

"""
    versions(job::Job)

Returs the valid versions of `job`.
"""
versions(job::Job) = job_versions(main_job_dir(job))
version(job::Job) = job.version
version(dir::AbstractString) = occursin(VERSION_DIR_NAME, dir) ? parse(Int, splitpath(dir)[end]) : main_job_version(dir)

function last_job_version(dir::AbstractString)
    versions = job_versions(dir)
    return isempty(versions) ? 0 : versions[end]
end

"""
    last_version(job::Job)

Returns the last version number of `job`.
"""
last_version(job::Job) = last_job_version(main_job_dir(job))

function version_dir(dir::AbstractString, version::Int)
    tpath = joinpath(dir, VERSION_DIR_NAME, "$version")
    return tpath
end
version_dir(job::Job) = version_dir(main_job_dir(job), job.version)
version_dir(job::Job, version::Int) = version_dir(main_job_dir(job), version)

exists_version(dir::AbstractString, version::Int) = version ∈ job_versions(dir)
exists_version(job::Job, version::Int) = exists_version(main_job_dir(job), version)

"""
    maybe_cp_main_version(job::Job)

Looks in the `job.dir` for the version of the job in the main directory, and copies it to the
respective directory in the `.versions`.
"""
function maybe_cp_main_version(job::Job)
    maindir = main_job_dir(job)
    if ispath(joinpath(maindir, "job.tt"))
        tjob = load_job(maindir)
        cp(tjob, joinpath(tjob, VERSION_DIR_NAME, "$(job.version)"); force = true)
    end
end

"""
    switch_version!(job::Job[, version::Int])

Switches the version of `job` to one of the previously stored ones.
It will save also the current version for future reference.
"""
function switch_version!(job::Job, version::Int)
    cur_version = job.version
    if version != cur_version
        version_assert(job, version)
        out = load_job(main_job_dir(job); version = version)
        for f in fieldnames(Job)
            if f == :server
                continue
            end
            setfield!(job, f, getfield(out, f))
        end
    end
    return job
end

function switch_version!(job::Job)
    vs = versions(job)
    timestamps = []
    for v in vs
        mdatapath = joinpath(version_dir(job, v), ".metadata.jld2")
        if ispath(mdatapath) &&
           haskey(load(mdatapath), "metadata") &&
           haskey(load(mdatapath)["metadata"], :timestamp)
            push!(timestamps,
                  string(round(load(mdatapath)["metadata"][:timestamp], Dates.Second)))
        else
            push!(timestamps, "")
        end
    end

    menu = RadioMenu(join.(zip(["v$v" for v in string.(vs)], timestamps),
                           ("\tsaved on: ",)))
    choice = request("Please select which version to switch to:", menu)
    if choice != -1
        return switch_version!(job, vs[choice])
    end
end

function version_assert(job, version)
    @assert exists_version(job, version) "Version $version does not exist for job."
end

"""
    rm_version!(job::Job, version::Int)
    rm_versions!(job::Job, versions::Int...)

Removes the specified `versions` from the `job` if they exist.
"""
function rm_version!(job::Job, version::Int)
    version_assert(job, version)
    if version == main_job_version(job)
        for f in readdir(job)
            if occursin(".workflow", f) || f == VERSION_DIR_NAME || splitext(f)[2] == ".jl"
                continue
            else
                rm(joinpath(job, f))
            end
        end
        md = main_job_dir(job)
        lv = last_version(job)
        if lv == version
            lv = versions(job)[end-1]
        end
        if lv != 0
            real_path = version_dir(md, lv)
            for f in readdir(real_path)
                cp(joinpath(real_path, f), joinpath(md, f), force=true)
            end
        end
    else
        rm(version_dir(job, version); recursive = true)
    end
end

"Removes temporary directories of the specified job versions."
function rm_tmp_dirs!(job, vers = versions(job)...)
    for v in vers
        version_assert(job, v)
        rm(joinpath(version_dir(job, v), TEMP_CALC_DIR); recursive = true)
    end
end


text:
using Documenter

# Print `@debug` statements (https://github.com/JuliaDocs/Documenter.jl/issues/955)
if haskey(ENV, "GITHUB_ACTIONS")
    ENV["JULIA_DEBUG"] = "Documenter"
end

using MCMCDiagnosticTools

DocMeta.setdocmeta!(
    MCMCDiagnosticTools, :DocTestSetup, :(using MCMCDiagnosticTools); recursive=true
)

makedocs(;
    modules=[MCMCDiagnosticTools],
    authors="David Widmann",
    repo="https://github.com/devmotion/MCMCDiagnosticTools.jl/blob/{commit}{path}#{line}",
    sitename="MCMCDiagnosticTools.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://devmotion.github.io/MCMCDiagnosticTools.jl",
        assets=String[],
    ),
    pages=["Home" => "index.md"],
    strict=true,
    checkdocs=:exports,
)

deploydocs(;
    repo="github.com/devmotion/MCMCDiagnosticTools.jl", push_preview=true, devbranch="main"
)


text:
# ---
# title: 1329. Sort the Matrix Diagonally
# id: problem1329
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Sort
# link: <https://leetcode.com/problems/sort-the-matrix-diagonally/description/>
# hidden: true
# ---
#
# Given a `m * n` matrix `mat` of integers, sort it diagonally in ascending
# order from the top-left to the bottom-right then return the sorted array.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/01/21/1482_example_1_2.png)
#
#
#
#     Input: mat = [[3,3,1,1],[2,2,1,2],[1,1,1,2]]
#     Output: [[1,1,1,1],[1,2,2,2],[1,2,3,3]]
#
#
#
#
# **Constraints:**
#
#   * `m == mat.length`
#   * `n == mat[i].length`
#   * `1 <= m, n <= 100`
#   * `1 <= mat[i][j] <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function large_neighborhood_search!(
    solution::Solution, mapf::MAPF; neighborhood_size=1, steps=10, progress=true
)
    @assert is_feasible(solution, mapf)
    cost = flowtime(solution, mapf)
    p = Progress(steps; desc="LNS steps: ", enabled=progress)
    for _ in 1:steps
        next!(p)
        agents = random_neighborhood(mapf, neighborhood_size)
        backup = remove_agents!(solution, agents)
        cooperative_astar!(solution, agents, mapf)
        new_cost = flowtime(solution, mapf)
        if is_feasible(solution, mapf) && new_cost < cost  # keep
            cost = new_cost
        else  # revert
            for a in agents
                solution[a] = backup[a]
            end
        end
    end
    return solution
end

function large_neighborhood_search(mapf::MAPF; neighborhood_size=1, steps=10, progress=true)
    solution = cooperative_astar(mapf, shuffle(1:nb_agents(mapf)))
    large_neighborhood_search!(
        solution, mapf; neighborhood_size=neighborhood_size, steps=steps, progress=progress
    )
    return solution
end


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

JLLWrappers.@generate_main_file_header("IOAPI")
JLLWrappers.@generate_main_file("IOAPI", UUID("21284ac2-4ede-594e-ad4f-4494079aca48"))
end  # module IOAPI_jll


text:
using PkgSRA, Test, Random
Random.seed!(13)

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

# Generate test data: Lorenz
include("../examples/example_lorenz.jl")
dat = Array(solve_lorenz_system())
numerical_grad = numerical_derivative(dat, ts)

# Generate chain
sindy_library = Dict("cross_terms"=>2,"constant"=>nothing);
my_sindy_model = sindyc(dat, numerical_grad,
                        library=sindy_library, use_lasso=true)
my_turing_model = convert_sindy_to_turing_enforce_zeros(my_sindy_model;
                                dat_noise_prior=Normal(0.0, 1.0),
                                coef_noise_std=0.01)

Turing.turnprogress(false)
chain = generate_chain(dat, numerical_grad, my_turing_model,
                            iterations=5,
                            num_training_pts=100, start_ind=101)[1]

# Number of terms tests
terms = get_nonzero_terms(my_sindy_model)
@test length(chain.name_map.parameters) - 1 == length(terms)
@test chain.name_map.parameters[end] == "noise"

# Converting back to SINDy model
my_sindy_sample = sindy_from_chain(my_sindy_model,
                                   chain,
                                   enforced_zeros=true)
terms_sample = get_nonzero_terms(my_sindy_sample)
@test terms == terms_sample

# Test helper prediction function
pred1 = my_sindy_model(dat[:,1], 0)
pred2 = sindy_predict(my_sindy_model,
                      my_sindy_model.A[terms],
                      dat[:,1],
                      nonzero_terms=terms)
@test all(pred1 .== pred2)


text:
# Define type
abstract type MeanSquaredDisplacement end

# Define structs
struct DisplacementPositions <: MeanSquaredDisplacement
     p::RandomWalker.Walker2D
     q::RandomWalker.Walker2D
end

struct TaxiCabDistance{T} <: MeanSquaredDisplacement
     d::T
end

struct TaxiCabDistanceXY{T} <: MeanSquaredDisplacement
    dₓ::T
    dy::T
end


text:
# ---
# title: 992. Subarrays with K Different Integers
# id: problem992
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Hash Table, Two Pointers, Sliding Window
# link: <https://leetcode.com/problems/subarrays-with-k-different-integers/description/>
# hidden: true
# ---
#
# Given an array `A` of positive integers, call a (contiguous, not necessarily
# distinct) subarray of `A` _good_ if the number of different integers in that
# subarray is exactly `K`.
#
# (For example, `[1,2,3,1,2]` has `3` different integers: `1`, `2`, and `3`.)
#
# Return the number of good subarrays of `A`.
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [1,2,1,2,3], K = 2
#     Output: 7
#     Explanation: Subarrays formed with exactly 2 different integers: [1,2], [2,1], [1,2], [2,3], [1,2,1], [2,1,2], [1,2,1,2].
#
#
# **Example 2:**
#
#
#
#     Input: A = [1,2,1,3,4], K = 3
#     Output: 3
#     Explanation: Subarrays formed with exactly 3 different integers: [1,2,1,3], [2,1,3], [1,3,4].
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 20000`
#   2. `1 <= A[i] <= A.length`
#   3. `1 <= K <= A.length`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
type Point2
  x::Int64
  y::Int64
  #z = 3
end

p = Point2(1, 2)
println(p.x)
p.x += 3
p.y += 2
println(p)


text:
module NanStats

# package code goes here
include("nan_statsjc.jl")

export nanmean, nanvar, nanskew, nankurt

end # module


text:
module CBF

import ..FileFormats
import MathOptInterface

const MOI = MathOptInterface

MOI.Utilities.@model(
    Model,
    (MOI.Integer,),
    (),
    (
        MOI.Reals,
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        MOI.DualExponentialCone,
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction,)
)

function MOI.supports_constraint(
    ::Model{T},
    ::Type{MOI.VariableIndex},
    ::Type{<:MOI.Utilities.SUPPORTED_VARIABLE_SCALAR_SETS{T}},
) where {T}
    return false
end

function MOI.supports_constraint(
    ::Model,
    ::Type{MOI.VariableIndex},
    ::Type{MOI.Integer},
)
    return true
end

function MOI.supports(
    ::Model,
    ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}},
)
    return false
end

"""
    Model()

Create an empty instance of `FileFormats.CBF.Model`.
"""
Model() = Model{Float64}()

Base.show(io::IO, ::Model) = print(io, "A Conic Benchmark Format (CBF) model")

include("read.jl")
include("write.jl")

end


text:
# push!(LOAD_PATH, "/home/ian/Grad/Research/Disks/DiskJockey/")
push!(LOAD_PATH, "/home/ian/Grad/Research/Disks/DiskJockey/movie/")
# push!(LOAD_PATH, "/n/home07/iczekala/DiskJockey/")
# push!(LOAD_PATH, "/pool/scout0/DiskJockey/")

using consonance

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--run_index", "-r"
    help = "Output run index"
    default = 0
    arg_type = Int
end

parsed_args = parse_args(ARGS, s)
run_index = parsed_args["run_index"]

imgdir = moviedir * @sprintf("g%02d/", run_index)
scratchdir = moviedir * "scratch/"

import YAML
config = YAML.load(open("config.yaml"))

# make the directory where RADMC will make all its files exist
if !ispath(imgdir)
    println("Creating ", imgdir)
    mkdir(imgdir)
end

if !ispath(scratchdir)
    mkdir(scratchdir)
end

cd(scratchdir)
println("Now in ", pwd())

files = ["lines.inp", "molecule_co.inp", "wavelength_micron.inp", "radmc3d.inp"]
for file in files
    src = moviedir * file
    dst = scratchdir * file
    cp(src, dst, remove_destination=true)
end

using DiskJockey.constants
using DiskJockey.image
using DiskJockey.model

species = "12CO"

function make_image(pars, id::Int)

    write_model(pars, "", grid, species)

    run(pipeline(`radmc3d image incl $(pars.incl) posang $(pars.PA) npix $npix loadlambda`, DevNull))

    src = "image.out"
    dst = imgdir * @sprintf("gimage%04d.out", id)

    cp(src, dst, remove_destination=true)

end

grd = config["grid"]

const global npix = config["npix"] # number of pixels
const global grid = Grid(grd["nr"], grd["ntheta"], grd["r_in"], grd["r_out"], true)
write_grid("", grid)
write_lambda(lams, "")

start = run_index * nframes_per_proc + 1
ids = Int[i for i=start:(start + nframes_per_proc)]

for i in ids
    if i <= nframes
        make_image(pars[i], i)
    else
        break
    end
end

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


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

struct AdminUsersSessionApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_adminUsersSessionInvalidate(_api::AdminUsersSessionApi, team_id::String, session_id::Int32; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.users.session.invalidate", ["slackAuth"])
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.form, "session_id", session_id)  # type Int32
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Invalidate a single session for a user by session_id
Param: team_id::String (required)
Param: session_id::Int32 (required)
Return: Dict{String, Any}
"""
function adminUsersSessionInvalidate(_api::AdminUsersSessionApi, team_id::String, session_id::Int32; _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersSessionInvalidate(_api, team_id, session_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminUsersSessionInvalidate(_api::AdminUsersSessionApi, response_stream::Channel, team_id::String, session_id::Int32; _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersSessionInvalidate(_api, team_id, session_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminUsersSessionReset(_api::AdminUsersSessionApi, user_id::String; mobile_only=nothing, web_only=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.users.session.reset", ["slackAuth"])
    Swagger.set_param(_ctx.form, "user_id", user_id)  # type String
    Swagger.set_param(_ctx.form, "mobile_only", mobile_only)  # type Bool
    Swagger.set_param(_ctx.form, "web_only", web_only)  # type Bool
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Wipes all valid sessions on all devices for a given user
Param: user_id::String (required)
Param: mobile_only::Bool
Param: web_only::Bool
Return: Dict{String, Any}
"""
function adminUsersSessionReset(_api::AdminUsersSessionApi, user_id::String; mobile_only=nothing, web_only=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersSessionReset(_api, user_id; mobile_only=mobile_only, web_only=web_only, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminUsersSessionReset(_api::AdminUsersSessionApi, response_stream::Channel, user_id::String; mobile_only=nothing, web_only=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersSessionReset(_api, user_id; mobile_only=mobile_only, web_only=web_only, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export adminUsersSessionInvalidate, adminUsersSessionReset


text:
using Distributed

# */ --------------------------------------------------------------------------------------------- /* #

function get_case()::Int64

    println("")
    println(" case 0 -> 0 uncertain parameters")
    println(" case 1 -> 1 uncertain parameters")
    println(" case 2 -> 2 uncertain parameters")
    println(" case 3 -> 3 uncertain parameters")
    println(" select case study:")
    cs = ""
    cs = readline()
    while (cs ∉ ["0","1","2","3"])
        println(" error: select 0, 1, 2, or 3")
        println(" select case:")
        cs = readline()
    end

    return parse(Int64,cs)
end

function get_algm()::Int64

    println("")
    println(" algorithm 0 -> deterministic equivalent")
    println(" algorithm 1 -> Stand_Bend (Standard Benders)")
    println(" algorithm 2 -> Adapt_Bend (Benders with Adaptive Oracles)")
    println(" algorithm 3 -> Zaker_Bend (Benders, Zakeri et al.)")
    println(" select Benders-type algorithm:")
    al = ""
    al = readline()
    while (al ∉ ["0","1","2","3"])
        println(" error: select 0, 1, 2, or 3")
        println(" select algorithm:")
        al = readline()
    end

    return parse(Int64,al)
end

function get_parw(cs::Int64)::Int64

    cs == 0 ? W = 2   : nothing
    cs == 1 ? W = 12  : nothing
    cs == 2 ? W = 90  : nothing
    cs == 3 ? W = 756 : nothing
    println("")
    println(" number of subproblems to solve each iter j:")
    println(" input w (0 < w < $W)")
    w = ""
    w = readline()
    while (parse(Float64,w) <= 0 || parse(Float64,w) >= W )
        println(" error: select 0 < w < $W")
        println(" input w:")
        w = readline()
    end

    return parse(Float64,w)
end

function get_parσ()::Float64

    println("")
    println(" select parameter σ:")
    println(" input σ (1 < σ <= 100)")
    σ = ""
    σ = readline()
    while (parse(Float64,σ) <= 1 || parse(Float64,σ) > 100 )
        println(" error: select 1 < σ <= 100")
        println(" input σ:")
        σ = readline()
    end

    return parse(Float64,σ)
end

function get_params()::Tuple{Int64,Int64,Float64,Int64}

    c = get_case()
    a = get_algm()
    a == 2 ? w = get_parw(c) : w =  1
    a == 3 ? q = get_parσ( ) : q = .0
    return c,a,q,w
end

# */ --------------------------------------------------------------------------------------------- /* #


text:
@testset "Results" begin
    # Inline definition of the Sphere function
    f_sphere(x) = sum(x.^2)

    RANDOM_SEED = 6918234
    rng = MersenneTwister(RANDOM_SEED)

    @testset "Benchmarks" begin
        total_iterations = 20000
        val = PSO(
            Easom(),
            Population(35, 2, -100.0, 100.0, rng),
            total_iterations,
            rng
        )
        ground_truth = @SVector [π, π]
        @test isapprox(val.x, ground_truth, atol = 1e-8)
        @test isapprox(val.min, -1.0)
        @test val.impl ≡ "PSO"
        @test val.iterations ≡ total_iterations
    end

    @testset "Functions" begin
        total_iterations = 20000
        val = PSO(
            f_sphere,
            Population(35, 3, -100.0, 100.0, rng),
            total_iterations,
            rng
        )
        @test isapprox(val.x, zeros(3), atol=eps())
        @test isapprox(val.min, 0.0, atol=eps())
        @test val.impl ≡ "PSO"
        @test val.iterations ≡ total_iterations
        optim_res = OptimizationResults(val.x,
                f_sphere(val.x),
                "PSO",
                val.iterations)
        @show optim_res
    end
end


text:
using CuArrays, CUDAnative, CUDAdrv

#This is the kernel. The plus and minus ones relative to the C version are due
#to the fact that ONE INDEXING SUCKS
function dot(a,b,c, N, threadsPerBlock, blocksPerGrid)

    #Set up shared memory cache for this current block. Has to be dynamic to be
    #able to take command line args - this will not work with
    #@cuStaticSharedMem
    cache = @cuDynamicSharedMem(Int64, threadsPerBlock)

    #Initialise some variables. Minus ones are to get the indexing consistent
    #with the C version
    tid = (threadIdx().x - 1) + (blockIdx().x - 1) * blockDim().x
    cacheIndex = threadIdx().x - 1
    temp::Int64 = 0

    #iterate over vector to do dot product in parallel way
    while tid < N
        temp += a[tid + 1] * b[tid + 1]
        tid += blockDim().x * gridDim().x
    end

    #set cache values
    cache[cacheIndex + 1] = temp

    #synchronise threads
    sync_threads()

    #In the step below, we add up all of the values stored in the cache
    i::Int = blockDim().x/2
    while i!=0
        if cacheIndex < i
            cache[cacheIndex + 1] += cache[cacheIndex + i + 1]
        end
        sync_threads()
        i/=2
    end

    #cache[1] now contains the sum of vector dot product calculations done in
    #this block, so we write it to c
    if cacheIndex == 0
        c[blockIdx().x] = cache[1]
    end

    return nothing
end

#Self explanatory
function sum_squares(x)
    return (x * (x + 1) * (2 * x + 1) / 6)
end

function main()

    #Can't type global variables in Julia, which messes things up.
    #So instead we locally define these and pass them explicitly as arguments to
    #the kernel
    N::Int64 = 33 * 1024
    threadsPerBlock::Int64 = 256
    blocksPerGrid::Int64 = min(32, (N + threadsPerBlock - 1) / threadsPerBlock)

    #Create a,b and c
    a = CuArrays.CuArray(fill(0, N))
    b = CuArrays.CuArray(fill(0, N))
    c = CuArrays.CuArray(fill(0, blocksPerGrid))

    #Fill and b
    for i in 1:N
        a[i] = i
        b[i] = 2*i
    end

    #Execute the kernel. Note the shmem argument - this is necessary to allocate
    #space for the cache we allocate on the gpu with @cuDynamicSharedMem
    @cuda blocks = blocksPerGrid threads = threadsPerBlock shmem =
    (threadsPerBlock * sizeof(Int64)) dot(a,b,c, N, threadsPerBlock, blocksPerGrid)

    #Copy c back from the gpu (device) to the host
    c = Array(c)

    local result = 0

    #Sum the values in c
    for i in 1:blocksPerGrid
        result += c[i]
    end

    #Check whether output is correct
    println("Does GPU value ", result, " = ", 2 * sum_squares(N - 1))
end

main()


text:
# needs a variable debugging defined to load

error("Outdated: Remove Rhat")

println("Time to load BITEModel.jl")
@time import BITEModel;
const BM=BITEModel
const I2=BM.ITMIX2
using VAWTools, KissMCMC
using FileIO, JLD
import DataStructures: OrderedDict

# update_cache = false
# verbose = false
# plotyes = false
# itmixversion = 2

const progress_bar = true
const save_results = true

"""
    setup_one_exp(exp, gl, gb, pp, pm, pn,
                      n_1d_btilde, n_1d_fsl, n_1d_temp)

Essentially calls `BM.init_inverse` but also sets up:
- theta0
- removes all unused radar-lines from gb

Also has the extra feature to make more experiments:
- exp0: no radar lines
- exp 17-32: no IV, otherwise as 1.16

Above two only make sense (in the sense that there is something to fit)
for glaciers with IV data.
"""
function setup_one_exp(exp, gl, gb, pp, pm, pn,
                       n_1d_btilde, n_1d_fsl, n_1d_temp,
                       runtyp, fit_target, rm_iv;
                       th0_with_sigma=false)
    gln = BM.ITMIX2.make_ITMIX2_glacier(gl, exp, rm_iv)
    gbn = BM.Bands(gb, gl=gln)

    th0d = BM.theta0_dict_defaults(pm, n_1d_btilde, n_1d_fsl, n_1d_temp,
                                   with_sigma=th0_with_sigma)

    # update IC for some glaciers with otherwise zero density of prior
    # because of terminus flux or continuous glacier prior.
    if gl.gid==BITEModel.ITMIXGlacier(:Mocho, 2)
        # Mocho has zero thickness at the top with the given btilde
        # (which leads to prior==-Inf)
        th0d[:btilde][1] = 1.5
    elseif gl.gid==BITEModel.ITMIXGlacier(:Austfonna, 2)
        th0d[:btilde][1] = 1.0
    elseif gl.gid==BITEModel.ITMIXGlacier(:Freya, 2)
        th0d[:btilde][1] = 0.5
    elseif gl.gid==BITEModel.ITMIXGlacier(:Kesselwandferner, 2)
        th0d[:btilde][1] = 0.5
    end
    return gln, gbn, th0d, BM.init_inverse(gbn, pp, pm, pn,
                                           runtyp=runtyp,
                                           theta0_dict=th0d,
                                           fit_target=fit_target
                                           )...
end

"""
    run_one_exp(exp, gl, gb, pp, pm, pn, pl;
                     n_1d_btilde = 5,
                     n_1d_fsl = 3,
                     n_1d_temp = 1)

Runs one ITMIX2 experiment.  Useful to use in `map, `pmap`, etc.

Also has the extra feature to make more experiments:
- exp0: no radar lines
- exp 17-32: no IV, otherwise as 1.16

Above two only make sense (in the sense that there is something to fit)
for glaciers with IV data.
"""
function run_one_exp(exp, gl, gb, pp, pm, pn, pl, runtyp, fit_target;
                     rm_iv=false,
                     n_1d_btilde=3,
                     n_1d_fsl=3,
                     n_1d_temp=1,
                     dir="results_Werder"
                     )
    # I get an odd bug in pmap with ERROR: LoadError: On worker 2: UndefVarError: h_iv not defined
    # when using symbols for runtyp and fit_target, thus the clutch of passing them as strings:
    runtyp = Symbol(runtyp)

    println("\n --- Running experiment $exp of $(BM.getname(gl))")

    (gln, gbn, th0d, theta0, logposterior, logposterior1d, logprior, logprior_debug,
     loglikelihood, loglikelihood1d, fwdm_fn, fwdm1d_fn,
     pmcmc_defaults, fit_target) =
         setup_one_exp(exp, gl, gb, pp, pm, pn,
                       n_1d_btilde, n_1d_fsl, n_1d_temp,
                       runtyp, fit_target, rm_iv)
    pmcmc = BM.MCMCNum(;pmcmc_defaults...)
    # some tests:
    if debugging
        fwdm_fn(theta0.th0);
        print("    Time to run forward model:")
        @time fwdsol = fwdm_fn(theta0.th0);
        print("    Posterior value:")
        @show logposterior(theta0.th0)[1]
        print("    Time to calculate posterior:")
        @time logposterior(theta0.th0)[1]
        print("    Starting MCMC:")
    end

    if logprior(theta0.th0)==-Inf
        logprior_debug(theta0.th0)
        error("Zero density of prior of $gl at $(theta0.th0)")
    end

    # sample posterior
    sol = BM.mcmc(logposterior, theta0, pmcmc;
                  verbose=debugging,
                  use_progress_meter=progress_bar,
                  misc=Dict(:logposterior=>logposterior,
                            :logposterior1d=>logposterior1d,
                            :logprior=>logprior,
                            :logprior_debug=>logprior_debug,
                            :loglikelihood=>loglikelihood,
                            :loglikelihood1d=>loglikelihood1d,
                            :fwdm_fn=>fwdm_fn,
                            :fit_target=>fit_target))
    debugging && println("MCMC of experiment $exp of $(BM.getname(gl)) done in $(signif(wall_time/60,4)) minutes")
    debugging && println("   Rhat=$(mean(sol.Rhat)), $(maximum(sol.Rhat))")
    #debugging && print_results(thetas, accept_ratio, names=BM.get_varnames(theta0))
    if save_results
        fln_jld = I2.save_itmix2_run(sol.mh2d, sol.sh2d, gln, pl, sol.miv2d, sol.siv2d, dir=dir)[end]
        FileIO.save(fln_jld,
                    "Rhat", sol.Rhat,
                    "sample_size", sol.sample_size,
                    "nthin", sol.nthin,
                    "thetas", sol.thetas,
                    "reshape_revert", sol.reshape_revert,
                    "runtyp", runtyp,
                    "fit_target", fit_target,
                    "theta0_names", sol.theta0.names,
                    "theta0_th0", sol.theta0.th0,
                    "th0_dict", th0d,
                    "wall_time", wall_time)
    end

    return res
end

println("itmix-setup.jl loaded")


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

nested_error_expr = quote
    try
        __not_a_binding__
    catch
        1 ÷ 0  # Generate error while handling error
    end
end

nested_error_pattern = r"""
    ERROR: DivideError: integer division error
    Stacktrace:.*

    caused by: UndefVarError: __not_a_binding__ not defined
    Stacktrace:.*
    """s

@testset "display_error" begin
    # Display of errors which cause more than one entry on the exception stack
    excs = try
        eval(nested_error_expr)
    catch
        Base.current_exceptions()
    end
    @test typeof.(first.(excs)) == [UndefVarError, DivideError]
    @test occursin(nested_error_pattern, sprint(Base.display_error, excs))

    @test occursin(r"""
        2-element ExceptionStack:
        DivideError: integer division error
        Stacktrace:.*

        caused by: UndefVarError: __not_a_binding__ not defined
        Stacktrace:.*
        """s, sprint(show, excs))
end

@testset "Fallback REPL" begin
    # Fallback REPL should show errors with display_error
    errio = IOBuffer()
    Base.eval_user_input(errio, nested_error_expr, true)
    err_str = String(take!(errio))
    @test occursin(nested_error_pattern, err_str)
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
using Knet
libknet8handle = Libdl.dlopen(Knet.libknet8)

SIZE1 = 1000
SIZE2 = 100
ITER = 100000
mat32b = rand(Float32,SIZE1,SIZE2)
mat32 = KnetArray(mat32b)
col32b = rand(Float32,SIZE1)
col32 = KnetArray(col32b)
row32b = rand(Float32,1,SIZE2)
row32 = KnetArray(row32b)
mat64b = rand(Float64,SIZE1,SIZE2)
mat64 = KnetArray(mat64b)
col64b = rand(Float64,SIZE1)
col64 = KnetArray(col64b)
row64b = rand(Float64,1,SIZE2)
row64 = KnetArray(row64b)

f32 = f64 = nothing

function cuda21test(fname, jname=fname, o...)
    println(fname)
    fcpu = eval(parse(jname))

    global f32 = Libdl.dlsym(libknet8handle, fname*"_32_21")
    print("gpu col32 ")
    @time cuda21rep(f32,mat32,col32)
    print("cpu col32 ")
    @time out32 = fcpu21rep(fcpu,mat32b,col32b)
    isapprox(out32,Array(col32)) || warn("$fname 32 mat col")
    print("gpu row32 ")
    @time cuda21rep(f32,mat32,row32)
    print("cpu row32 ")
    @time out32 = fcpu21rep(fcpu,mat32b,row32b)
    isapprox(out32,Array(row32)) || warn("$fname 32 mat row")

    global f64 = Libdl.dlsym(libknet8handle, fname*"_64_21")
    print("gpu col64 ")
    @time cuda21rep(f64,mat64,col64)
    print("cpu col64 ")
    @time out64 = fcpu21rep(fcpu,mat64b,col64b)
    isapprox(out64,Array(col64)) || warn("$fname 64 mat col")
    print("gpu row64 ")
    @time cuda21rep(f64,mat64,row64)
    print("cpu row64 ")
    @time out64 = fcpu21rep(fcpu,mat64b,row64b)
    isapprox(out64,Array(row64)) || warn("$fname 64 mat row")
end

function cuda21rep{T}(f,x::KnetArray{T},y::KnetArray{T})
    nx = length(x); ny = length(y)
    if size(y,1)==1; sy=size(x,1); else; sy=1; end
    for i=1:ITER
        ccall(f,Void,(Cint,Ptr{T},Cint,Cint,Ptr{T}),nx,x,sy,ny,y)
    end
    Knet.@cuda(cudart,cudaDeviceSynchronize,())
    Knet.@cuda(cudart,cudaGetLastError,())
end

function fcpu21rep{T}(f,x::Array{T},y::Array{T})
    if size(y,1)==1; region=1; else; region=2; end
    out = nothing
    for i=1:1 # ITER
        out = f(x,region)
    end
    return out
end

for f in Knet.cuda21
    isa(f,Tuple) || (f=(f,))
    cuda21test(f...)
end


text:
# Contains all the types definitions

using Distributions

abstract type SobolEvals end

struct SobolDists
    P::Int
    dist::Array{Distributions.Distribution,1}
    function SobolDists(P,dist)
        length(dist) != P ?
        error("Dimensions are incorrect!") :
        new(P,dist)
    end
end

mutable struct SobolSamples
    N::Int
    P::Int
    dist::Array{Distributions.Distribution,1}
    A::Array{Float64,2}
    B::Array{Float64,2}
    function SobolSamples(N,P,dist,A,B)
        (size(A) != (N,P) || size(B) != (N,P) || length(dist) != P) ?
        error("Dimensions are incorrect!") :
        new(N,P,dist,A,B)
    end
end

mutable struct SobolEvalsSingle <: SobolEvals
    N::Int
    P::Int
    dist::Array{Distributions.Distribution,1}
    fA::Array{Float64,2}
    fB::Array{Float64,2}
    fC::Array{Array{Float64,2}}
    function SobolEvalsSingle(N,P,dist,fA,fB,fC)
        (size(fA,1) != N || size(fB,1) != N || size(fC,1) != P) ?
        error("Function evaluation numbers don't match!") :
        new(N,P,dist,fA,fB,fC)
    end
end

mutable struct SobolEvalsJoint <: SobolEvals
    N::Int
    P::Int
    dist::Array{Distributions.Distribution,1}
    Cpairs::Array{Array{Int,1}}
    fC::Array{Array{Float64,2}}
    function SobolEvalsJoint(N,P,dist,Cpairs,fC)
        (size(fC,1) != binomial(P,2)) ?
        error("Function evaluation numbers don't match!") :
        new(N,P,dist,Cpairs,fC)
    end
end

mutable struct SobolEvalsSingleAndJoint <: SobolEvals
    N::Int
    P::Int
    dist::Array{Distributions.Distribution,1}
    fA::Array{Float64,2}
    fB::Array{Float64,2}
    fC::Array{Array{Float64,2}}
    Cpairs::Array{Array{Int,1}}
    fCij::Array{Array{Float64,2}}
    function SobolEvalsSingleAndJoint(N,P,dist,fA,fB,fC,Cpairs,fCij)
        (size(fA,1) != N || size(fB,1) != N || size(fC,1) != P || size(fCij,1) != binomial(P,2)) ?
        error("Function evaluation numbers don't match!") :
        new(N,P,dist,fA,fB,fC,Cpairs,fCij)
    end
end

max_stars_repo_name: mstobb/VarianceBasedSensitivity