text:
function optimize(
    objective,
    method;
    constraints = nothing,
    population = nothing,
    options::Options = Options(; default_options(method)...),
    state = nothing,
)::OptimizationResult

    if constraints === nothing
        constraints = NoConstraints()
    end

    if isa(objective, Function)
        objective = Objective(objective, first(population))
    end

    if population === nothing
        population = initial_population(method, objective)
    end

    if state === nothing
        state = initial_state(method, objective, population, options)
    end

    iteration = 0
    stopped = false
    converged, counter_tol = false, 0 # tolerance convergence
    is_moo = ismultiobjective(objective)
    start_time = stop_time = time()

    trace = []
    # Ignore first state
    # if options.store_trace
    #     push!(trace, copy(state))
    # end

    # 终止条件
    # 1. 收敛
    # 2. 达到指定的 iteration
    # 3. 达到指定的运行时间
    # 4. terminate(state)?
    # 5. callback

    while !converged && !stopped && iteration < options.iterations
        iteration += 1
        state.iteration = iteration

        state.start_time = time()

        should_break = update_state!(
            method,
            state,
            objective,
            constraints,
            options
        )

        state.stop_time = time()

        converged = assess_convergence!(state)
        counter_tol = converged ? counter_tol + 1 : 0
        converged = converged && (counter_tol > options.successive_f_tol)
        converged = converged || terminate(state)

        if (options.store_trace)
            push!(trace, copy(state))
        end
        if (options.state_callback !== nothing)
            options.state_callback(state)
        end

        should_break && break
    end

    stop_time = time()

    return OptimizationResult(
        method,
        minimizer(state),
        is_moo ? NaN : value(state),
        iteration,
        iteration == options.iterations,
        converged,
        state.metrics,
        f_calls(objective),
        trace,
        options.time_limit,
        stop_time - start_time,
        is_moo,
    )
end


text:
function pointstoforces(
    edges::Array{GeometricalPredicates.Line2D{GeometricalPredicates.Point2D},1},
    scaler::Scaler,
    fscale::Float64,
    /pfix,
    fh::Function,
)::Dict{Point2D,Array{Float64,1}}
    bars = Array{Point2D,1}()
    barvec = Array{Array{Float64,1},1}()
    pointstofvec = Dict{Point2D,Array{Float64,1}}()
    for edge in edges
        b = unscaledpoint2d(getb(edge), scaler)
        a = unscaledpoint2d(geta(edge), scaler)
        push!(bars, barcenter(a, b))
        push!(barvec, barvector(a, b))
        push!(pointstofvec, geta(edge) => [0, 0])
        push!(pointstofvec, getb(edge) => [0, 0])
    end
    fvec = forcevectors(bars, barvec, fscale, fh)
    return pointstoforces(edges, pointstofvec, fvec, pfix)
end

function barcenter(a::Point2D, b::Point2D)::Point2D
    return Point(getx(a) + ((getx(b) - getx(a)) / 2), gety(a) + ((gety(b) - gety(a)) / 2))
end

function barvector(a::Point2D, b::Point2D)::Array{Float64, 1}
    return [getx(b) - getx(a), gety(b) - gety(a)]
end

function forcevectors(
    bars::Array{Point2D,1},
    barvec::Array{Array{Float64,1},1},
    fscale::Float64,
    fh::Function,
)
    L = [sqrt(sum(vsum .^ 2)) for vsum in barvec]
    hbars = [fh(getx(p), gety(p)) for p in bars]
    L0 = hbars * fscale * sqrt(sum(L .^ 2) / sum(hbars .^ 2))
    fvec = maximum.(L0 - L) ./ L .* barvec
end

function pointstoforces(
    edges::Array{GeometricalPredicates.Line2D{GeometricalPredicates.Point2D},1},
    pointstofvec::Dict{Point2D,Array{Float64,1}},
    fvec,
    pfix
)::Dict{Point2D,Array{Float64,1}}
    iterator = 1
    for edge in edges
        a = geta(edge)
        b = getb(edge)
        push!(pointstofvec, a => pointstofvec[a] + (-fvec[iterator]))
        push!(pointstofvec, b => pointstofvec[b] + (fvec[iterator]))
        iterator += 1
    end
    for p in eachrow(pfix)
        p = Point(p[1], p[2])
        delete!(pointstofvec, p)
        push!(pointstofvec, p => [0, 0])
    end
    return pointstofvec
end

ext: jl

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
#=
1. Construct N parallel cuboids of size 10000 × 1 × 100 spaced one unit apart in y -direction as object W .
2. Construct N parallel cuboids of size 1 × 10000 × 100 spaced one unit apart in x -direction as object W ′ .
3. Align W and W ′ at their lower front left corner.
4. Move W ′ along the z-axis for fifty units.
5. Measure time for W ∪ W′.
=#

using ViewerGL; GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation;

n = 30
V1,(_,EV1,FV1,_) = Lar.cuboid([1,2n,10],true);
cuboid1 = (V1,FV1,EV1);
V2,(_,EV2,FV2,_) = Lar.cuboid([2n,1,10],true);
cuboid2 = (V2,FV2,EV2);

# @time let
	group1 = Lar.struct2lar(Lar.Struct(repeat([cuboid1, Lar.t(2,0,0)],outer=n)));
	group2 = Lar.struct2lar(Lar.Struct(repeat([cuboid2, Lar.t(0,2,0)],outer=n)));
	assembly = Lar.Struct([group1, Lar.t(0,0,5), group2]);

	V,FV,EV = Lar.struct2lar(assembly);
	cop_EV = convert(Lar.ChainOp, Lar.coboundary_0(EV::Lar.Cells));
	cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells); ## TODO: debug
	W = convert(Lar.Points, V');
	V, copEV, copFE, copCF = Lar.space_arrangement( W, cop_EV, cop_FE);
	W = convert(Lar.Points, V');
	V,CVs,FVs,EVs = Lar.pols2tria(W, copEV, copFE, copCF);
# end;

GL.VIEW( GL.GLExplode(V,CVs[2:end],5,5,5,99,0.5) );
GL.VIEW( GL.GLExplode(V,CVs[1:1],5,5,5,99,0.5) );


text:
# new Sandshark example
# add more julia processes
# using Distributed
# nprocs() < 4 ? addprocs(4-nprocs()) : nothing

using Caesar, RoME
# @everywhere using Caesar, RoME
using Interpolations
using Distributions

using RoMEPlotting
using Gadfly, DataFrames
using ProgressMeter
using DelimitedFiles
using Logging

# const TU = TransformUtils

Gadfly.set_default_plot_size(35cm,25cm)

include(joinpath(@__DIR__,"Plotting.jl"))
include(joinpath(@__DIR__,"SandsharkUtils.jl"))

# Step: Selecting a subset for processing and build up a cache of the factors.
epochs = timestamps[50:2:60]
lastepoch = 0
for ep in epochs
  global lastepoch
  if lastepoch != 0
    # @show interp_yaw(ep)
    deltaAng = interp_yaw(ep) - interp_yaw(lastepoch)

    wXi = TU.SE2([interp_x(lastepoch);interp_y(lastepoch);interp_yaw(lastepoch)])
    wXj = TU.SE2([interp_x(ep);interp_y(ep);interp_yaw(ep)])
    iDXj = se2vee(wXi\wXj)
    NAV[ep] = iDXj
    # println("$(iDXj[1]), $(iDXj[2]), $(iDXj[3])")

    odoDict[ep] = Pose2Pose2(MvNormal(NAV[ep], Matrix(Diagonal([0.1;0.1;0.005].^2))))
  end
  rangepts = rangedata[ep][:]
  rangeprob = kde!(rangepts)
  azipts = azidata[ep][:,1]
  aziprob = kde!(azipts)

  # prep the factor functions
  ppbrDict[ep] = Pose2Point2BearingRange(aziprob, rangeprob)
  lastepoch = ep
end

## Step: Building the factor graph
fg = initfg()
# Add a central beacon with a prior
addVariable!(fg, :l1, Point2)
# Pinger location is (0.6; -16)
addFactor!(fg, [:l1], PriorPose2( MvNormal([0.6; -16], Matrix(Diagonal([0.1; 0.1].^2)) ) ), autoinit=false)

index = 0
for ep in epochs
    global index
    curvar = Symbol("x$index")
    addVariable!(fg, curvar, Pose2)
    if ep != epochs[1]
      # Odo factor x(i-1) -> xi
      addFactor!(fg, [Symbol("x$(index-1)"); curvar], odoDict[ep], autoinit=false)
    else
      # Prior to the first pose location (a "GPS" prior)
      initLoc = [interp_x(ep);interp_y(ep);interp_yaw(ep)]
      println("Adding a prior at $curvar, $initLoc")
      addFactor!(fg, [curvar], PriorPose2( MvNormal(initLoc, Matrix(Diagonal([0.1;0.1;0.05].^2))) ), autoinit=false)
    end
    index+=1
end

# Just adding the first one...
# addFactor!(fg, [:x0; :l1], ppbrDict[epochs[1]], autoinit=false)
addFactor!(fg, [:x5; :l1], ppbrDict[epochs[6]], autoinit=false)

# Debugging options
getSolverParams(fg).drawtree = true
getSolverParams(fg).showtree = true
getSolverParams(fg).async = true
getSolverParams(fg).downsolve = false
getSolverParams(fg).multiproc = false
getSolverParams(fg).limititers = 30

tree = solveTree!(fg, recordcliqs=ls(fg))

drawPosesLandms(fg, meanmax=:max)

0

# assignTreeHistory!(tree, hist)
#
# csmAnimate(fg, tree, ls(fg), frames=1000) #[:x0; :x2]
# # Base.rm("/tmp/caesar/csmCompound/out.ogv")
# run(`ffmpeg -r 10 -i /tmp/caesar/csmCompound/csm_%d.png -c:v libtheora -vf fps=25 -pix_fmt yuv420p -vf "scale=trunc(iw/2)*2:trunc(ih/2)*2" -q 10 /tmp/caesar/csmCompound/out.ogv`)
# @async run(`totem /tmp/caesar/csmCompound/out.ogv`)

# getVariableInferredDim(fg, :x0)
# getVariableInferredDim(fg, :x1)
# getVariableInferredDim(fg, :x2)
# getVariableInferredDim(fg, :x3)
# getVariableInferredDim(fg, :x4)
# getVariableInferredDim(fg, :x5)
# getVariableInferredDim(fg, :l1)

## OR Do it manually??
tree = wipeBuildNewTree!(fg)
drawTree(tree, show=true)

smt, hist = solveCliq!(fg, tree, :x0)
smt, hist = solveCliq!(fg, tree, :l1)

smt, hist = solveCliq!(fg, tree, :x4)
smt, hist = solveCliq!(fg, tree, :x2)
smt, hist = solveCliq!(fg, tree, :x1)

### DEV work below

cliq = getClique(tree,:x0)

# getCliqVariableInferredPercent(fg, cliq)

### Pretend solve cliq 2

cliq = getClique(tree, :x2)
syms = getCliqAllVarSyms(fg, cliq)
c3sfg = buildSubgraphFromLabels(fg, syms)
# c2sfg = buildSubgraphFromLabels(fg, syms)

dfg = fg
@info "8a, needs down message -- attempt down init"
prnt = getParent(tree, cliq)[1]

# take atomic lock when waiting for down ward information
# lockUpStatus!(getData(prnt))

dwinmsgs = prepCliqInitMsgsDown!(dfg, tree, prnt, logger=ConsoleLogger()) # cliqSubFg
dwnkeys = collect(keys(dwinmsgs))

## DEVIdea
msgfcts = addMsgFactors!(c3sfg, dwinmsgs)
# writeGraphPdf(c3sfg, show=true)

@info "8a, attemptCliqInitD., dwinmsgs=$(dwnkeys)"

# DEVidea
# sdims = getCliqVariableMoreInitDims(c3sfg, cliq) # removed after IIF v0.25
updateCliqSolvableDims!(cliq, sdims)

# determine if more info is needed for partial
# priorize solve order for mustinitdown with lowest dependency first
# follow example from issue #344
mustwait = false
if length(intersect(dwnkeys, getCliqSeparatorVarIds(cliq))) == 0 # length(dwinmsgs) == 0 ||
  @info "8a, attemptCliqInitDown_StateMachine, no can do, must wait for siblings to update parent first."
  global mustwait = true
elseif getSiblingsDelayOrder(tree, cliq, prnt, dwinmsgs, logger=ConsoleLogger())
  @info "8a, attemptCliqInitD., prioritize"
  global mustwait = true
elseif getCliqSiblingsPartialNeeds(tree, cliq, prnt, dwinmsgs, logger=ConsoleLogger())
  @info "8a, attemptCliqInitD., partialneedsmore"
  global mustwait = true
end

# remove the downward messages too
deleteMsgFactors!(c3sfg, msgfcts)

# unlock
# @info "8a, attemptCliqInitD., unlockUpStatus!"
# unlockUpStatus!(getData(prnt))

fetchCliqSolvableDims(getClique(tree, :x4))
fetchCliqSolvableDims(getClique(tree, :x2))
fetchCliqSolvableDims(getClique(tree, :x0))

prnt = getClique(tree, :x1)
getCliqSiblingsPriorityInitOrder(tree, prnt)

# getData(getClique(tree, :x0)).solvableDims

# fetch(getData(getClique(tree, :x0)).solvableDims)

#


text:
using BikeStats
using Test

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


text:
using MPI
using HDF5
using Test

@testset "mpio" begin

MPI.Init()

info = MPI.Info()
comm = MPI.COMM_WORLD

nprocs = MPI.Comm_size(comm)
myrank = MPI.Comm_rank(comm)

@test HDF5.has_parallel()

let fileprop = p_create(HDF5.H5P_FILE_ACCESS)
    HDF5.h5p_set_fapl_mpio(fileprop, comm, info)
    h5comm, h5info = HDF5.h5p_get_fapl_mpio(fileprop)

    # check that the two communicators point to the same group
    if isdefined(MPI, :Comm_compare)  # requires recent MPI.jl version
        @test MPI.Comm_compare(comm, h5comm) === MPI.CONGRUENT
    end
end

# open file in parallel and write dataset
fn = MPI.bcast(tempname(), 0, comm)
A = [myrank + i for i = 1:10]
h5open(fn, "w", comm, info) do f
    @test isopen(f)
    g = g_create(f, "mygroup")
    dset = d_create(g, "B", datatype(Int64), dataspace(10, nprocs), chunk=(10, 1), dxpl_mpio=HDF5.H5FD_MPIO_COLLECTIVE)
    dset[:, myrank + 1] = A
end

MPI.Barrier(comm)
h5open(fn, comm) do f  # default: opened in read mode, with default MPI.Info()
    @test isopen(f)
    @test keys(f) == ["mygroup"]

    B = read(f, "mygroup/B", dxpl_mpio=HDF5.H5FD_MPIO_COLLECTIVE)
    @test !isempty(B)
    @test A == vec(B[:, myrank + 1])

    B = f["mygroup/B", dxpl_mpio=HDF5.H5FD_MPIO_COLLECTIVE]
    @test !isempty(B)
    @test A == vec(B[:, myrank + 1])
end

MPI.Barrier(comm)

B = h5read(fn, "mygroup/B", fapl_mpio=(comm, info), dxpl_mpio=HDF5.H5FD_MPIO_COLLECTIVE)
@test A == vec(B[:, myrank + 1])

MPI.Barrier(comm)

B = h5read(fn, "mygroup/B", (:, myrank + 1), fapl_mpio=(comm, info), dxpl_mpio=HDF5.H5FD_MPIO_COLLECTIVE)
@test A == vec(B)

# we need to close HDF5 and finalize the info object before finalizing MPI
finalize(info)
HDF5.h5_close()

MPI.Barrier(MPI.COMM_WORLD)

MPI.Finalize()

end # testset mpio


text:
module autograd

include("./variable.jl")
include("./helper.jl")

import .variable: Variable, Node
import .helper: ensure_float

function backward(v::Variable)

    function compute_grads(v::Variable, child_grad::Float64)

        for node in v.parent_nodes

            parent_var = node.variable
            if parent_var.requires_grad
                parent_grad = node.grad * child_grad
                parent_var.grad += parent_grad
            else
                continue
            end

        compute_grads(parent_var, parent_grad)
        end

    end
    compute_grads(v, 1.)

end

function zero_grad(v::Variable)
    v.grad = 0.
end

export backward, zero_grad

end # Module

max_stars_repo_name: brymer-meneses/Autograd

max_issues_repo_name: brymer-meneses/Autograd

max_forks_repo_name: brymer-meneses/Autograd


text:
struct Model
    fw
    fc
    fs
end

function Model(wordembeds::Matrix{T}, charembeds::Matrix{T}, ntags::Int) where T
    fw = @graph n begin
        Node(Lookup(wordembeds), n)
    end

    d = size(charembeds, 1)
    fc = @graph n begin
        n = Node(Lookup(charembeds), n)
        n = Node(Conv1D(T,5d,5d,2d,d), n)
        Node(max, n, 2)
    end

    d = size(wordembeds,1) + size(charembeds,1)*5
    fs = @graph n begin
        n = Node(BiLSTM(T,d,d), n)
        Node(Linear(T,2d,ntags), n)
    end
    Model(fw, fc, fs)
end

function (m::Model)(word::Var, chars::Vector{Var})
    w = m.fw(word)
    cs = Var[]
    for i = 1:length(chars)
        push!(cs, m.fc(chars[i]))
    end
    c = cat(2, cs...)
    s = cat(1, w, c)
    m.fs(s)
end

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
#=
  Bayesian Data Analysis of Iris dataset

  This model uses a linear regression.

  Compared to iris.jl this program splits the data into a training and test set.
  Here are some results for different training percentages.

  * For train_pct 0.75
    Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

           a    1.0930    0.2342     0.0074    0.0087   383.1404    1.0002       90.8346
   params[1]   -0.0735    0.0656     0.0021    0.0029   352.8822    1.0033       83.6610
   params[2]   -0.0659    0.0661     0.0021    0.0031   404.0677    1.0012       95.7960
   params[3]    0.2289    0.0625     0.0020    0.0026   303.2456    1.0009       71.8932
   params[4]    0.5658    0.1042     0.0033    0.0036   374.0015    0.9993       88.6680
       sigma    0.2228    0.0155     0.0005    0.0005   630.9561    0.9995      149.5866

    Num correct: 37.0
    Num incorrect: 1.0
    Accuracy: 0.9736842105263158
    Confusion matrix:
    3×3 Matrix{Float64}:
    15.0   0.0  0.0
     0.0  14.0  0.0
     0.0   1.0  8.0
    num bad: 1
    Distribution of the training categories:
    1.0 : 35 (31.25%)
    2.0 : 36 (32.142857142857146%)
    3.0 : 41 (36.607142857142854%)

    * train pct 0.66
    Num correct: 49.0
    Num incorrect: 2.0
    Accuracy: 0.9607843137254902
    Confusion matrix:
    3×3 Matrix{Float64}:
    18.0   0.0   0.0
    0.0  15.0   1.0
    0.0   1.0  16.0
    num bad: 2
    Distribution of the training categories:
    1.0 : 32 (32.323232323232325%)
    2.0 : 34 (34.343434343434346%)
    3.0 : 33 (33.333333333333336%)

    * train pct 0.066 (10 train instances)
    Num correct: 130.0
    Num incorrect: 10.0
    Accuracy: 0.9285714285714286
    Confusion matrix:
    3×3 Matrix{Float64}:
    46.0   0.0   0.0
     0.0  42.0   6.0
     0.0   4.0  42.0
    num bad: 10
    Distribution of the training categories:
    1.0 : 4 (40.0%)
    2.0 : 2 (20.0%)
    3.0 : 4 (40.0%)

    Not too bad.

    * train pct 0.02 (3 traingin instances)
    Interestingly, even with as few as 3 (of 150) training instances the
    model might get quite good accuracy, but sometimes it's really bad,
    e.g. when the training instances are of the same class.
    Here are some runs with some different distributions of the three
    instances:

    Num correct: 135.0
    Num incorrect: 12.0
    Accuracy: 0.9183673469387755
    Confusion matrix:
    3×3 Matrix{Float64}:
    49.0   0.0   0.0
     0.0  38.0  11.0
     0.0   1.0  48.0
    num bad: 12
    Distribution of the training categories:
    1.0 : 1 (33.333333333333336%)
    2.0 : 1 (33.333333333333336%)
    3.0 : 1 (33.333333333333336%)

    Num correct: 112.0
    Num incorrect: 35.0
    Accuracy: 0.7619047619047619
    Confusion matrix:
    3×3 Matrix{Float64}:
    49.0   0.0   0.0
     0.0  48.0   0.0
     0.0  35.0  15.0
    num bad: 35
    Distribution of the training categories:
    1.0 : 1 (33.333333333333336%)
    2.0 : 2 (66.66666666666667%)

    Num correct: 121.0
    Num incorrect: 26.0
    Accuracy: 0.8231292517006803
    Confusion matrix:
    3×3 Matrix{Float64}:
    48.0   0.0   0.0
     0.0  24.0  26.0
     0.0   0.0  49.0
    num bad: 26
    Distribution of the training categories:
    1.0 : 2 (66.66666666666667%)
    3.0 : 1 (33.333333333333336%)

    Num correct: 144.0
    Num incorrect: 3.0
    Accuracy: 0.9795918367346939
    Confusion matrix:
    3×3 Matrix{Float64}:
    50.0   0.0   0.0
     0.0  46.0   3.0
     0.0   0.0  48.0
    num bad: 3
    Distribution of the training categories:
    2.0 : 1 (33.333333333333336%)
    3.0 : 2 (66.66666666666667%)

    Num correct: 55.0
    Num incorrect: 92.0
    Accuracy: 0.3741496598639456
    Confusion matrix:
    3×3 Matrix{Float64}:
    47.0   0.0  0.0
    42.0   8.0  0.0
    20.0  30.0  0.0
    num bad: 92
    Distribution of the training categories:
    1.0 : 3 (100.0%)

    Cf iris.jl
=#
using Turing, StatsPlots, Distributions, StatsBase
using DataFrames,RDatasets
include("jl_utils.jl")

# Read dataset
iris = dataset("datasets", "iris")

num_rows, num_cols = size(iris)

# train_pct = 0.75
# train_pct = 0.66
# train_pct = 0.066 # train with 10 instances of 150. Not too bad result!
train_pct = 0.02 # train with 3 instances. Sometimes quite good results.

# Convert Categorical to integer
conv = Dict(
    "setosa" => 1,
    "versicolor" => 2,
    "virginica" => 3,
)
display(describe(iris))
println()

iris2 = Matrix(transform(iris,[:Species] => ByRow(i->conv[i]))[:,[1,2,3,4,6]])
x = iris2[:,1:4]
y = iris2[:,5]

println("\nMean value for each species: ", keys(conv))
display(keys(conv) .|> s -> (s,mean(x[iris.Species.==s,:], dims=1)))
println("\n")

x_train, x_test, y_train, y_test = split_train_test(x,y,train_pct)
num_train_rows = length(y_train)
num_test_rows = length(y_test)

println("Num train instances: $num_train_rows")
println("Num test instances: $num_test_rows")

println("Distribution of the training categories:")
for (k,v) in sort(collect(make_hash(y_train)))
  println("$k : $v ($(100*v/num_train_rows)%)")
end

model = linear_regression_model(x_train,y_train)

# println("Prior:")
# chns = sample(model, Prior(), 10_000)
# display(chns)

println("\nModel:")
# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 1_000)
# chns = sample(model, PG(5), 1_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)
# chns = sample(model, HMC(0.1,6), 1_000)
chns = sample(model, NUTS(), 1_000)

display(chns)

println("\nparameter: credible_intervals")
for v in ["a", "params[1]", "params[2]", "params[3]", "params[4]","sigma"]
  credible_interval(chns,v,0.95)
end

# display(plot(chns))
println("\nPredictions:")

mpred = linear_regression_model(x_test,missing)
pred = predict(mpred,chns)
display(pred)

cats = Dict(1=>1, 2=>2, 3=>3)
confusion, bads = confusion_matrix(y_test,pred,cats,false)
println("Confusion matrix:")
display(confusion)
# println("\nBad classfied:")
# display(bads)
println("num bad: $(length(bads))")
#=
println("\nbads: credible_intervals")
for b in bads
  println("b:",b)
  ix = b[1]
  println("x[$ix]: $(x[ix,:])  true cat (y_test[$ix]): $(y_test[ix])  (pred mean: ", mean(pred["y[$ix]"]))
  credible_interval(pred,"y[$(b[1])]",0.3)
end
=#

println("Distribution of the training categories:")
for (k,v) in sort(collect(make_hash(y_train)))
  println("$k : $v ($(100*v/num_train_rows)%)")
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
Mach2_benchmark = OneDOneQMachine(0.0, #R
                        1.3125, #Xd
                        1.2578, #Xq
                        0.1813, #Xd_p
                        0.25, #Xq_p
                        5.89, #Td0_p
                        0.6, #Tq0_p
                        100.0)   #MVABase

Mach3_benchmark = OneDOneQMachine(0.0, #R
                        1.3125, #Xd
                        1.2578, #Xq
                        0.1813, #Xd_p
                        0.25, #Xq_p
                        5.89, #Td0_p
                        0.6, #Tq0_p
                        100.0)   #MVABase

Shaft2_benchmark = SingleMass(3.01, #H (M = 6.02 -> H = M/2)
                       0.0) #D

Shaft3_benchmark = SingleMass(3.01, #H (M = 6.02 -> H = M/2)
                      0.0) #D

no_pss = PSSFixed(0.0)

fixed_tg = TGFixed(1.0) #eff

gen2_avr_benchmark_m = AVRTypeI(20.0, #Ka - Gain
                   0.01, #Ke
                   0.063, #Kf
                   0.2, #Ta
                   0.314, #Te
                   0.35, #Tf
                   0.001, #Tr
                   5.0, #Vrmax
                   -5.0, #Vrmin
                   0.0039, #Ae - 1st ceiling coefficient
                   1.555) #Be - 2nd ceiling coefficient

gen3_avr_benchmark_m = AVRTypeI(20.0, #Ka - Gain
                                0.01, #Ke
                                0.063, #Kf
                                0.2, #Ta
                                0.314, #Te
                                0.35, #Tf
                                0.001, #Tr
                                5.0, #Vrmax
                                -5.0, #Vrmin
                                0.0039, #Ae - 1st ceiling coefficient
                                1.555) #Be - 2nd ceiling coefficient

Gen2_benchmark = DynGenerator(1, #number
                              :Gen1, #name
                              nodes_benchmark[2],
                              1.0, #ω_ref
                              1.0142, #V_ref
                              1.0, #P_ref
                              Mach2_benchmark, #machine
                              Shaft2_benchmark, #shaft
                              gen2_avr_benchmark_m, #avr
                              fixed_tg, #tg
                              no_pss) #pss

Gen3_benchmark = DynGenerator(2, #number
                                :Gen2, #name
                                nodes_benchmark[3],
                                1.0, #ω_ref
                                1.0059, #V_ref
                                1.0, #P_ref
                                Mach3_benchmark, #machine
                                Shaft3_benchmark, #shaft
                                gen3_avr_benchmark_m, #avr
                                fixed_tg, #tg
                                no_pss) #pss

inf_gen_benchmark = StaticSource(1, #number
                 :InfBus, #name
                 nodes_benchmark[1], #bus
                 1.02, #VR
                 0.0, #VI
                 0.000001) #Xth

ThreeBus_Benchmark = DynamicSystem(nodes_benchmark, #buses
                                  branch_benchmark, #branches
                                  [Gen2_benchmark, Gen3_benchmark], #dynamic injections
                                  [inf_gen_benchmark, loads_benchmark[1], loads_benchmark[2], loads_benchmark[3]],
                                  100.0, #Sbase
                                  60.0) #fbase


text:
using BenchmarkTools
using LuxurySparse
using LinearAlgebra
using SparseArrays

id = IMatrix(1<<8)
v = randn(1<<8) + im*randn(1<<8)
Dv = Diagonal(v)

#sp = sprand(ComplexF64, 1<<8, 1<<8, 0.03)
#pm = PermuteMultiply(randperm(1<<8), randn(1<<8))
sp = SparseMatrixCSC(Dv)
ds = rand(ComplexF64, 1<<4,1<<4)
pm = PermMatrix(Dv)

mats = [id, sp, ds, pm, Dv]
for A in mats
    bg = BenchmarkGroup()
    for B in mats
        if !(A===B===ds)
            println("============ $(typeof(A))-$(typeof(B)) =============")
            display(@benchmark kron($A, $B))
            println()
        end
    end
    display(run(bg, verbose=true))
end


text:
"""
    taxon(df::DataFrame, id::Integer)

Returns a fully formed `NCBITaxon` based on its id. The `name` of the taxon
will be the valid scientic name associated to this id.
"""
function taxon(df::DataFrame, id::Integer)
    matched_indices = findall(isequal(id), df.tax_id)
    isempty(matched_indices) && throw(IDNotFoundInBackbone(id))
    submatches = df[matched_indices, :]
    @assert NCBITaxonomy.class_scientific_name in submatches.class
    sciname_position = findfirst(
        isequal(NCBITaxonomy.class_scientific_name), submatches.class
    )
    return NCBITaxon(submatches.name[sciname_position], id)
end

"""
    taxon(id::Integer)

Performs a search in the entire taxonomy backbone based on a known ID.
"""
taxon(id::Integer) = taxon(NCBITaxonomy.names_table, id)

function _id_from_name(name::AbstractString; kwargs...)
    return _id_from_name(NCBITaxonomy.names_table, name; kwargs...)
end

function _id_from_name(
    df::DataFrame,
    name::AbstractString;
    strict::Bool=true,
    dist::Type{SD}=Levenshtein,
    casesensitive::Bool=true,
    rank::Union{Nothing,Symbol}=nothing,
) where {SD<:StringDistance}
    if !isnothing(rank)
        @assert rank ∈ unique(df.rank)
        df = df[findall(isequal(rank), df.rank),:]
    end
    if strict
        positions = if casesensitive
            findall(isequal(name), df.name)
        else
            findall(isequal(lowercase(name)), df.lowercase)
        end
        # If the array is empty, we throw the "no name" error
        isempty(positions) && throw(NameHasNoDirectMatch(name))
        # If the array has a single element, this is the ticket
        length(positions) == 1 && return df.tax_id[first(positions)]
        # If neither of these are satisfied, the name has multiple matches and we throw the appropriate error
        taxa = taxon.(df.tax_id[positions])
        throw(NameHasMultipleMatches(name, taxa))
    else
        correct_name, position = findnearest(name, df.name, dist())
        return df.tax_id[position]
    end
end

"""
    taxon(name::AbstractString; kwargs...)

The `taxon` function is the core entry point in the NCBI taxonomy. It takes a
string, and a series of keywords, and go look for this taxon in the dataframe
(by default the entire names table).

The keywords are:

- `strict` (def. `true`), allows fuzzy matching
- `dist` (def. `Levenshtein`), the string distance function to use
- `casesensitive` (def. `true`), whether to strict match on lowercased names
- `rank` (def. `nothing`), the taxonomic rank to limit the search
"""
taxon(name::AbstractString; kwargs...) = taxon(NCBITaxonomy.names_table, name; kwargs...)

"""
    taxon(df::DataFrame, name::AbstractString; kwargs...)

Additional method for `taxon` with an extra dataframe argument, used most often
with a `namefinder`. Accepts the usual `taxon` keyword arguments.
"""
function taxon(df::DataFrame, name::AbstractString; kwargs...)
    id = _id_from_name(df, name; kwargs...)
    isnothing(id) && return nothing
    return taxon(df, id)
end

"""
    ncbi_str(s)

A string macro to perform a strict taxonomic search.
"""
macro ncbi_str(s)
    return taxon(s)
end


text:
function iterative_barycenter_initialization(measurements::Vector{Pointcloud}
                                            ,weights::Vector{R}
                                            ,c::S
                                            ,eps::Float64
                                            ) where {R <: Real, S <: Real}
    #ospa_iterative_sampling_space(measurements, weights, c, eps)
    ospa_barycenter_sol2(measurements, weights, c, eps)
end

function ospa_barycenter_sol2(measurements::Vector{Pointcloud},
                              weights::Vector{R},
                              c::S,
                              eps::Float64
                             ) where {R <: Real, S <: Real}
    #assignments = Vector(map(x -> zeros(Integer, size(x)[1]), measurements))
    assignments = Vector(map(x -> fill(0, size(x)[1]), measurements))
    first_guess = init_by_point(measurements, assignments)
    ospa_barycenter_sol2(measurements, assignments, weights, c, first_guess, eps)
end

#function ospa_barycenter_sol2(measurements::Vector{Pointcloud},
#                              assignments::Vector{Vector{Integer}},
#                              weights::Vector{R},
#                              c::S,
#                              init_first_guess::Function,
#                              ;
#                              do_postprocessing=false,
#                              confidence_threshold=0.05,
#                             ) where {R <: Real, S <: Real}
#    # Start with first guess
#    ospa_barycenter_sol2(measurements, assignments, weights, c, init_fun, bc)
#end

function ospa_iterative_sampling_space(measurements::Vector{Pointcloud}
                                       ,weights::Vector{R}
                                       ,c::S
                                       ,eps::Float64
                                      ) where {R <: Real, S <: Real}
    cost_old = Inf
    #eps = 0.05
    assignments = []
    converged = false
    while !converged
        bc, cost, assignments, subcluster_cost = update_barycenter_sol2(bc, measurements, weights, c)
        if cost_old - cost < eps
            converged = true
        end
        cost_old = cost
    end

    bc, cost, assignments, subcluster_cost = add_remaining_points(bc, measurements, assignments, weights, c, subcluster_cost)
    sorting = sortperm(subcluster_cost)
    bc = bc[sorting, :]
    subcluster_cost = subcluster_cost[sorting, :]
    subcluster_cost = subcluster_cost / sum(subcluster_cost)
    bc, subcluster_cost
end

function ospa_barycenter_sol2(measurements::Vector{Pointcloud},
                              assignments::Vector{Vector{I}},
                              weights::Vector{R},
                              c::S,
                              bc::Pointcloud,
                              eps::Float64
                             ) where {R <: Real, S <: Real, I <: Integer}
    cost_old = Inf
    #eps = 0.05
    assignments = []
    converged = false
    while !converged
        bc, cost, assignments, subcluster_cost = update_barycenter_sol2(bc, measurements, weights, c)
        if cost_old - cost < eps
            converged = true
        end
        cost_old = cost
    end

    # Old
    #bc, cost, assignments, subcluster_cost = update_barycenter_sol2(bc, measurements, weights, c)
    #println("first bc")
    #println(bc)
    bc, cost, assignments, subcluster_cost = add_remaining_points(bc, measurements, assignments, weights, c, subcluster_cost)
    #if do_postprocessing
    #    println(subcluster_cost)
    #    bc, cost, subcluster_cost = confidence_based_postprocessing(bc, measurements, assignments, cost, confidence_threshold, c, subcluster_cost)
    #end
    #now select optimal barycenter.
    #println("Subcluster cost: $(subcluster_cost)")
    sorting = sortperm(subcluster_cost)
    bc = bc[sorting, :]
    subcluster_cost = subcluster_cost[sorting, :]
    cost_rising = false
    first_ind = 1
    last_ind = 1
    previous_cost = Inf
    current_cost = Inf
    best_bc = bc[first_ind:last_ind, :]
    best_cost = mospa(best_bc, measurements, weights, c)
    #println("Now select best barycenter.")
    while !cost_rising && last_ind <= size(bc)[1]
        current_cost = mospa(bc[first_ind:last_ind, :], measurements, weights, c)
#        println("Current_cost: $(current_cost).")
 #       println(bc[first_ind:last_ind, :])
        #println("Current cost for size $(last_ind): $(current_cost).\n Previous cost: $(previous_cost).")
        if current_cost > previous_cost || last_ind >= size(bc)[1]
        #if last_ind >= size(bc)[1]
            cost_rising = true
        else
            last_ind += 1
            previous_cost = current_cost
        end
    end
    best_bc = bc[first_ind:last_ind, :]
    best_bc, current_cost
end

function add_remaining_points(bc::Pointcloud,
                              measurements::Vector{Pointcloud},
                              assignments::Vector{Array{Integer,1}},
                              weights::Vector{R},
                              c::S,
                              subcluster_cost::Vector{Float64},
                             ) where {R <: Real, S <: Real}
    if map(x -> any(x .== 0), assignments) |> any
        exhausted = false
    else
        exhausted = true
    end
    idx_next_sub_cluster = size(bc)[1]
    loop_counter = 1
    while !exhausted
        #println("loop_counter: $(loop_counter).")
        loop_counter += 1
        idx_next_sub_cluster += 1
        next_point, outer_ind, inner_ind = draw_from_remaining_pool(measurements, assignments)
        #println("Next point drawn: $(next_point), outer_ind: $(outer_ind), inner_ind: $(inner_ind).")
        #next_sub_cluster = zeros(next_point)
        #old_next_sub_cluster = zeros(next_point)
        next_sub_cluster = fill(0, size(next_point))
        old_next_sub_cluster = fill(0, size(next_point))
        tmp_assignment_log = []
        tmp_sub_cluster_cost = 0
        converged = false
        while !converged
            opt_ass = selective_optimal_assignments(next_point, measurements, assignments)
            #println("assignment: $(opt_ass).")
            new_weight = 0
            tmp_assignment_log = []
            next_sub_cluster = fill(0, size(next_point))
            tmp_sub_cluster_cost = 0
            for i=1:length(opt_ass)
                idx = opt_ass[i][1]
                if idx != 0
                    #println("Next point: $(next_point), meas: $(measurements[i][idx:idx, :]), dist: $(p2dist(next_point, measurements[i][idx:idx, :])[1]).")
                    if p2dist(next_point, measurements[i][idx:idx, :])[1] <= c
                        next_sub_cluster += measurements[i][idx:idx, :]*weights[i]/max(length(measurements), length(assignments[i]))
                        new_weight += weights[i]/max(length(measurements), length(assignments[i]))
                        push!(tmp_assignment_log, (i, idx, idx_next_sub_cluster))
                        tmp_sub_cluster_cost += 1
                    end
                end
            end
            if new_weight > 0
                next_sub_cluster /= new_weight
            end
            if p2dist(next_sub_cluster, old_next_sub_cluster)[1] == 0
                converged = true
            else
                old_next_sub_cluster = next_sub_cluster
            end
        end
        bc = vcat(bc, next_sub_cluster)
        #println("Next bc:")
        #println(bc)
        push!(subcluster_cost, tmp_sub_cluster_cost)
        #println("New assignments: $(tmp_assignment_log).")
        copy_new_assignments!(assignments, tmp_assignment_log)
        #println("Size: $(size(bc)).")
        exhausted = check_for_remaining_points(assignments)
    end
    bc, mospa(bc, measurements, weights, c), assignments, subcluster_cost
end

function confidence_based_postprocessing(bc::Pointcloud,
                                         measurements::Vector{Pointcloud},
                                         assignments::Vector{Vector{Integer}},
                                         cost::R,
                                         confidence_threshold::Real,
                                         c::S,
                                         subcluster_cost::Vector{Float64},
                                         ;
                                         niter=20
                                        ) where {S <: Real, R <: Real}
    #Maximum sized barcyenter computed. Calc the confidence parameters.
    confidences, points_in_radius = compute_confidence(bc, measurements, c)
    clusters_to_shuffle = find(confidences .< confidence_threshold)
    # Find all clusters that need to be reshuffled now.
    clusters_to_shuffle_final = copy(clusters_to_shuffle)
    for i=1:length(points_in_radius)
        for j=1:length(points_in_radius[i])
            if any(k in points_in_radius[i][j] for k in clusters_to_shuffle)
                clusters_to_shuffle_final = unique([points_in_radius[i][j]; clusters_to_shuffle_final])
            end
        end
    end
    #unassign all points that need to be reshuffled.
    for cluster in clusters_to_shuffle_final
        subcluster_cost[cluster] = 0
    end
    for i =1:length(assignments)
        for j=1:length(assignments[i])
            if assignments[i][j] in clusters_to_shuffle_final
                assignments[i][j] = 0
            end
        end
    end
    #Now recompute the barycenter
    bc_best = copy(bc)
    cost_best = copy(cost)
    for _=1:niter
        bc_new, cost_new, _, subcluster_cost = add_remaining_points(bc, measurements, assignments, weights, c, subcluster_cost)
        if cost_new < cost_best
            bc_best = bc_new
            cost_best = cost_new
        end
    end
    bc_best, cost_best, subcluster_cost
end

function compute_confidence(bc::Pointcloud,
                            measurements::Vector{Pointcloud},
                            c::S
                           ) where {S <: Real}
    #Compute for eacht point to how many sub clusters it can belong
    counters =  map(x -> fill(Vector{Integer}(), size(x)[1]), measurements)
    #point_confidences = [ones(size(counters[i])) for i=1:length(counters)]
    #point_counter_for_clusters = [zeros(size(counters[i])) for i=1:length(counters)]
    cluster_confidence = zeros(size(bc)[1])
    subclusters_cardinality = zeros(size(bc)[1])
    for i=1:length(measurements)
        dmat = p2dist(measurements[i], bc)
        for j=1:size(measurements[i])[1]
            points_in_radius = find(dmat[j,:] .< c)
            counters[i][j] = [counters[i][j] ; points_in_radius]
            for k=1:length(points_in_radius)
                cluster_confidence[points_in_radius[k]] += 1/length(counters[i][j])
                #point_counter_for_clusters[points_in_radius[k]] += 1
                subclusters_cardinality[points_in_radius[k]] += 1
            end
        end
    end
    #println(point_counter_for_clusters)
    #println(cluster_confidence)
    for i=1:length(cluster_confidence)
        #cluster_confidence[i] /= point_counter_for_clusters[i]
        cluster_confidence[i] /= subclusters_cardinality[i]
    end
    cluster_confidence, counters
end

function copy_new_assignments!(assignments::Vector{Vector{Integer}}, new_assignments)
    for triple in new_assignments
        assignments[triple[1]][triple[2]] = triple[3]
    end
end

function update_barycenter_sol2(cur_bc::Pointcloud,
                                measurements::Vector{Pointcloud},
                                weights::Vector{R},
                                c::S,
                               ) where {R <: Real, S <: Real}
    #updated_bc = Matrix{eltype(cur_bc)}(0,size(cur_bc)[2])
    updated_bc = Matrix{eltype(cur_bc)}(fill(0, 0, size(cur_bc)[2]))
    subcluster_weight = Vector{Float64}()
    assignments = optimal_assignments(cur_bc, measurements)
    assigned_clusters = map(x -> zeros(Integer, size(x)[1]), measurements)
    #not sure if weights are handled correctly here
    for idxp = 1:size(cur_bc)[1]
        assigned_points = zeros(eltype(cur_bc), 1,size(cur_bc)[2])
        new_weight = 0
        tmp_subcluster_weight = 0
        for i=1:length(assignments)
            idxa = assignments[i][idxp]
            if idxa != 0 && p2dist(cur_bc[idxp:idxp,:], measurements[i][idxa:idxa,:])[1] <= c
                # Is the correct weight index idxa or i???
                tmp_subcluster_weight += 1
                #is it correct to change that line?
                #assigned_points += measurements[i][idxa:idxa,:]*weights[idxa]/max(length(measurements), length(assignments[i]))
                assigned_points += measurements[i][idxa:idxa,:]*weights[i]/max(length(measurements), length(assignments[i]))
                new_weight += weights[i]/max(length(measurements), length(assignments[i]))
                assigned_clusters[i][idxa] = idxp
            end
        end
        if new_weight > 0
            push!(subcluster_weight, tmp_subcluster_weight)
            updated_bc = vcat(updated_bc, assigned_points /new_weight)
        end
    end
    if size(updated_bc)[1] == 0
        println(cur_bc)
        error("update_barycenter: This should not happen!")
    end
    updated_bc, mospa(updated_bc, measurements, weights, c), assigned_clusters, subcluster_weight
end

function init_by_measurement(measurements::Vector{Pointcloud},
                             rest...
                            )
    measurements[rand(1:length(measurements))]
end

function init_by_point(measurements::Vector{Pointcloud},
                       assignments::Vector{Vector{I}},
                      ) where {I <: Integer}
    outer_ind = map(x -> any(x .== 0), assignments) |> findall |> rand
    inner_ind = (assignments[outer_ind] .== 0) |> findall |> rand
    measurements[outer_ind][inner_ind:inner_ind, :]
end

function draw_random_points(base_set, lower, upper, rand_scaling)
    n = rand(lower:upper)
    idxs = randperm(size(base_set)[1])
    sample = base_set[idxs[1:n], :]
    sample + rand(size(sample))*rand_scaling
end

function draw_from_remaining_pool(measurements::Vector{Pointcloud}, assignments::Vector{Vector{Integer}})
    outer_ind = map(x-> any(x.==0), assignments) |> findall |> rand
    inner_ind = (assignments[outer_ind] .== 0) |> findall |> rand
    measurements[outer_ind][inner_ind:inner_ind, :], outer_ind, inner_ind
end

function selective_optimal_assignments(pc::Pointcloud, measurements::Vector{Pointcloud}, assignments::Vector{Vector{Integer}})
    opt_assignments = []
    for i=1:length(measurements)
        # println("Index: $(i), pc: $(pc), meas: $(measurements[i]), ass: $(assignments[i]).")
        dmat = restricted_p2dist(pc, measurements[i], assignments[i])
        if all(dmat .== Inf) || isempty(dmat) # The last part is ok?
            push!(opt_assignments, [0])
        else
            push!(opt_assignments, hungarian(dmat)[1])
        end
    end
    opt_assignments
end

function restricted_p2dist(x::Pointcloud, y::Pointcloud, assignments::Vector{Integer})
    dmat = p2dist(x,y)
    dmat[:, assignments .!= 0] .= Inf
    dmat
end

function check_for_remaining_points(not_assigned_points::Vector{Vector{Integer}})
    map(x -> all(x .!= 0), not_assigned_points) |> all
end


text:
function main(args)
    return args
end
# Load inline arguments for worker unit
num_workers, this = main(ARGS)
num_workers = parse(Int64, num_workers)
this        = parse(Int64, this)

#
# Important directories
simdir          = "/remote/ceph2/group/gerda/data/mpik/gerda-simulations/gerda-gems-db/pss-dev"
simdir_old      = "/remote/ceph/group/gerda/data/simulation/gerda-mage-sim"
mapping_file    = simdir_old*"/UTILS/det-data/ged-mapping.json"
parameters_file = simdir_old*"/UTILS/det-data/ged-parameters.json"
config_dir      = "../../../2020-02-06_8380701d_st_ffischer/res/Impurity_Scan/config-dep/"

#
# Load packages and functions
include("../src/init.jl")
include("../src/fct.jl")
include("../src/worker_fct.jl")

#
# Output paths and filter settings
plots_base_path = "../../waveforms/sim/plots/raw/"
base_path_raw   = "../../waveforms/sim/raw_wf/"
base_path       = "../../waveforms/sim/raw_wf/"
sampling_time   = 1u"ns";
no_new_files    = 0
parameters = JSON.parsefile("../dicts/electronics_parameters.json")

channels = []
for c in 0:1:36
    files = glob(base_path_raw * lpad(c, 2, "0") * "-" * channel_to_name[c] * "/cl-wf/raw*/*.h5");
    if length(files) > 0
        push!(channels, c)
    end
end
channels = get_share_for_worker(channels, num_workers, this)

for ch in channels
    Base.run(`clear`)
    @info("Start Ch" * lpad(ch, 2, "0") * " | " * channel_to_name[ch])
    @info(">---------------------------------------------------------------<")

    files = glob(base_path_raw * lpad(ch, 2, "0") * "-" * channel_to_name[ch] * "/cl-wf/raw*/*.h5")
    pro = Progress(length(files), dt=0.5,
        barglyphs=BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
        barlen=10)
    for file in files
        file_nr = findall(x->x==file, files)[1]
        filename  = split(file, basename(file))[1]
        filename *= channel_to_name[ch] * "_w_filter" * split(basename(file), channel_to_name[ch])[2]
        filename = split(filename, "/cl-wf/")[1] * "/w_filter/" * split(filename, "/cl-wf/")[2]
        if !isfile(filename)
            data_raw = HDF5.h5open(file, "r") do h5f
                LegendHDF5IO.readdata(h5f, "generated_waveforms")
            end
            first_half = 1:Int(size(data_raw,1) / 2)
            waveform = add_baseline_and_extend_tail.(data_raw.waveform[first_half], 6000, 12000)

            #
            # Reduce time resolution to fit data
            waveform_w_bl = []
            map(x-> push!(waveform_w_bl, RDWaveform(x.time[1:10:end], x.value[1:10:end])), waveform)
            GBP = parameters[channel_to_name[ch]]["par"]["GBP"]*1e6
            tau = parameters[channel_to_name[ch]]["par"]["tau"]*1e-6
            Cd  = parameters[channel_to_name[ch]]["par"]["Cd"]*1e-12
            Cf  = parameters[channel_to_name[ch]]["par"]["Cf"]*1e-12

            data = Table(
                energy       = [],
                multiplicity = [],
                timestamp    = [],
                run          = [],
                channel      = [],
                waveform     = []
            )
            for i in eachindex(waveform_w_bl)
                pulse = waveform_w_bl[i].value
                if sum(pulse) > 0
                    filtered_pulse = applyElectronics(pulse; Ts = 10e-9, GBP = GBP, tau = tau, Kv = 150e3, Cd = Cd, Cf = Cf, Rf = 500e6)
                    inter = find_intersect(filtered_pulse, maximum(filtered_pulse)/2, 5)
                    cut   = filtered_pulse[(inter-399):1:(inter+400)]
                    time  = (0:10:7990)u"ns"
                    cut_pulse = RDWaveform(time, cut)
                    append!(data, Table(
                            energy       = [data_raw[i].edep],
                            multiplicity = [data_raw[i].event_mult],
                            timestamp    = [data_raw[i].evtno],
                            run          = [0],
                            channel      = [unique(data_raw[i].detno)[1]],
                            waveform     = [cut_pulse]
                        )
                    )
                end

            end
            events = Table(
                energy       = VectorOfArrays(Array{typeof(data[1].energy),1}(data.energy)),
                multiplicity = Int32.(data.multiplicity),
                timestamp    = Int32.(data.timestamp),
                run          = Int32.(data.run),
                channel      = Int32.(data.channel),
                waveform     = ArrayOfRDWaveforms( Array{typeof(data[1].waveform), 1}(data.waveform) )
            )

            !isdir(dirname(filename)) ? mkpath(dirname(filename)) : ""
            HDF5.h5open(filename, "w") do h5f
                LegendHDF5IO.writedata(h5f, "data", events)
            end
        end
        next!(pro)
    end
end


text:
"""
	  simplexintersection(S1::Array{Float64, 2}, S2::Array{Float64, 2};
        tolerance::Float64 = 1/10^10) -> Float64

Computes the volume of intersection between two n-dimensional simplices
by boundary triangulation. The simplices `S1` and `S2` are arrays of
(n, n+1), where each column is a vertex.

Note: the returned volume is not corrected. It should be divided by factorial(dim)
to obtain the true volume.

## How are intersections computed?
Intersections are computed as follows:

1. Find minimal set of points generating the intersection volume. These points form
a convex polytope Pᵢ.
2. Triangulate the faces of Pᵢ into simplices.
3. Combine each boundary simplex with an interior point in Pᵢ. The set of
all such combinations form a triangulation of Pᵢ.
4. Calculate the volume of each simplex in the resulting triangulation. The
sum of these volumes is the volume of the intersection.
"""
function simplexintersection(simplex1, simplex2; tol::Float64 = 1/10^10)
  # Nasty hack until we re-do everything for vectors of vectors.
  S1 = Array(simplex1)
  S2 = Array(simplex2)

  # Dimension
  n = size(S1, 1)

  # Centroid and radii
  c1 = Circumsphere(S1)[2:n+1]
  c2 = Circumsphere(S2)[2:n+1]
  r1 = Circumsphere(S1)[1]
  r2 = Circumsphere(S2)[1]

  # Orientation of simplices
  orientation_S1 = det([ones(1, n + 1); S1])
  orientation_S2 = det([ones(1, n + 1); S2])

  if abs(orientation_S1) < tol || abs(orientation_S2) < tol
    return 0.0
  end

  # Set volume to zero initially. Change only if there is intersection
  IntVol = 0.0

# -------------------------------------
# Simplices intersect in some way
# -------------------------------------

  # If the (distance between centroids)^2-(sum of radii)^2 < 0,
  # then the simplices intersect in some way.
  dist_difference::Float64 = (transpose(c1 - c2) * (c1 - c2) - (r1 + r2)^2)[1]

  if dist_difference < 0
    # Find the number of points of each simplex contained within the
    # circumsphere of the other simplex

    vertices1InCircum2 = SomeVertexInCircumsphere(S1, r2, c2)
    vertices2InCircum1 = SomeVertexInCircumsphere(S2, r1, c1)

    # At least one circumsphere contains vertices of the other simplex
    #println("")

    if vertices1InCircum2 + vertices2InCircum1 >= 1
      βs1in2, βs2in1, ordered_vertices1, ordered_vertices2, numof1in2, numof2in1 =
        BarycentricCoordinates(S1,S2,orientation_S1,orientation_S2,tol)
      # Trivial intersections
      TriviallyContained = heaviside0([numof1in2 numof2in1] .- (n+1))
      IsSomeContained = sum(TriviallyContained, dims=2)[1]

      if IsSomeContained == 2.0 # The simplices coincide
        IntVol = abs(orientation_S1)
      elseif IsSomeContained == 1.0 # One simplex is contained in the other

        if TriviallyContained[1] == 1.0 # Simplex1 is contained in Simplex2
          IntVol = abs(orientation_S1)
        else # Simplex2 is contained in Simplex1
          IntVol = abs(orientation_S2)
        end
      else # No simplex contains the other

        #print("SharedVertices\t\t\t")

        Ncomm, ordered_vertices1, ordered_vertices2 = SharedVertices(βs1in2,ordered_vertices1,ordered_vertices2,numof1in2,numof2in1)

        # Is there any shared face?
        if Ncomm == n
          IntVol = SharedFaceVolume(S2, βs1in2, ordered_vertices1, ordered_vertices2)
        else # The simplices do not share a face.
          IntVert, ConvexExpIntVert = IntersectionOfBoundaries_NoStorage(S1,S2,βs1in2,βs2in1, ordered_vertices1, ordered_vertices2, numof1in2, numof2in1, Ncomm, tol)
          if !isempty(IntVert)
            IntVert,ConvexExpIntVert = PolytopeGeneratingVertices(S1,S2,IntVert,ConvexExpIntVert,βs1in2,βs2in1,ordered_vertices1,ordered_vertices2,numof1in2,numof2in1,Ncomm);
            IntVol = VolumeComputation(IntVert, ConvexExpIntVert)
          end
        end
      end
    else
      #println("No circumsphere of either simplex contains vertices of the other simplex")
    end
  end

  return IntVol
end

export simplexintersection


text:
include(joinpath(@__DIR__, "../vb-cb/plot_vb_results.jl"))  # To precompile

using Distributed
addprocs(4)

@everywhere include(joinpath(@__DIR__, "../vb-cb/plot_vb_results.jl"))

@everywhere function parse_elbos(f; nlast=100)
  ms = match.(r"(?<=elbo: )-\d+\.\d+", f)
  elbos = map(m -> parse(Float64, m.match), filter(m -> !isnothing(m), ms))
  length(elbos) > nlast || (nlast = 1)
  return mean(elbos[(end-nlast):end])
end

@everywhere function separatesamples(Y)
  sample_ids = sort(unique(Y.sample_id))
  return Matrix.(Y[Y.sample_id .== i, Not(:sample_id)]
                 for i in sample_ids)
end

@everywhere function get_attr(path)
  p = match(r"(?<=p)\d+\.\d+", path).match
  seed = match(r"(?<=seed)\d+", path).match
  return (p=parse(Float64, p),
          seed=parse(Int, seed))
end

@everywhere function postprocess(path)
  # Read CB data
  y, markernames = let
    p, seed = get_attr(path)
    data_path = "../data/cb_transformed_reduced_p$(p).csv"
    Y = coalesce.(CSV.read(data_path), NaN)
    (separatesamples(Y), names(Y))
  end

  println("Processing: $(path)")
  makeplots(y=y, path_to_output=path,
            markernames=markernames, nlam=30)
end

function get_best_seed_for_p(paths, p; nlast=100)
  p_paths = filter(path -> match(Regex(string(p)), path) != nothing, paths)
  elbos = [let
             f = open(readlines, joinpath(dirname(path), "log.txt"))
             seed = parse(Int, match(r"(?<=seed)\d+", path).match)
             elbo = parse_elbos(f, nlast=nlast)
             (elbo, seed)
           end for path in p_paths]
  best_idx = argmax([elbo[1] for elbo in elbos])
  best_seed = elbos[best_idx][2]
  return best_seed
end

### MAIN ###
# NOTE: Modify this.
# RESULTS_DIR = "results/vb-cb-psens-bs500"
RESULTS_DIR = "results/vb-cb-psens-bs2000"

# List of paths to results.
paths = [joinpath(root, file)
         for (root, dirs, files) in walkdir(RESULTS_DIR)
         for file in files
         if file == "output.bson"]

# ps used.
ps = unique([parse(Float64, m.match) for m in match.(r"(?<=p)\d+\.\d+", paths)])
best_seeds = Dict(p => get_best_seed_for_p(paths, p) for p in ps)
best_paths = ["$(RESULTS_DIR)/p$(p)_seed$(s)/output.bson"
              for (p, s) in best_seeds]

# Post process in parallel
@time status = pmap(postprocess, best_paths, on_error=identity)
status_sanitized = map(s -> s == nothing ? "success" : s, status)

# Printing success / failure status of runs.
for (path, s) in zip(best_paths, status_sanitized)
  println("$(path) => status: $(s)")
end

println("ALL DONE!")

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
const FW_VERSION=v"0.3.0-alpha"
#=
+-----x
|5 1 6
|4   2
|8 3 7
y
=#
const direx=[0,1,0,-1,-1,1,1,-1],
direy=[-1,0,1,0,-1,-1,1,1],
dirl=[4,1,2,3,8,5,6,7],
dirr=[2,3,4,1,6,7,8,5],
diroppo=[3,4,1,2,7,8,5,6]

const Float=Float64,
Target=Union{Entity,Block},
DrawContext=GraphicsContext


text:
using PyPlot

# sample problem
  n = [64, 128, 256, 512, 1024, 2048, 4096];
cpu = [63.376e-3, 103.985e-3, 292.813e-3, 872.886e-3, 3.957, 22.645, 79.344 ];
gpu = [18.923e-3, 18.869e-3, 24.897e-3, 71.766e-3, 238.528e-3, 925.671e-3,4.984   ];

# FourierFlows TwoDTurb problem
  n = [128, 256, 512, 1024, 2048, 4096]
cpu = [0.446, 1.130, 3.344, 19.654, 98.487, 376.201]
gpu = [0.104, 0.152, 0.414, 1.405, 5.459, 23.724]

fig, axs = subplots(ncols=1, nrows=2, figsize=(6, 7))

subplot(2,1,1)
loglog(n, cpu, "*:", label="CPU", linewidth=0.8,  markersize=8)
loglog(n, gpu, "o:", label="GPU", linewidth=0.8,  markersize=6)
loglog(n, 0.0001*n.^2, "-k", label="slope 2", linewidth=1)
ylabel("walltime for 100 time-steps [s]", fontsize=12)
title("TwoDTurb module benchmark", fontsize=16)
# xlim(70, 10000)
legend()

subplot(2,1,2)
semilogx(n, cpu./gpu, "d:k", linewidth=0.8, markersize=6)
ylabel("GPU over CPU speedup", fontsize=12)
xlabel("number of grid-points per spatial dimension", fontsize=12)
# xlim(20, 10000)
ylim(0, 20)
legend()

savename = string("CPUversusGPU")
savefig(savename, dpi=300)


text:
include("Experiment.jl")

using JuLIP, BenchmarkTools
using Experiment

A = [1.0 sin(π/6); 0.0 cos(π/6)]
# n = Int64(10.0 / 1.0 + 1.0)
n = 50
t = 1:n
o = ones(n)
xref = A * [(t * o')[:]'; (o * t')[:]']
xvecs = [xref; zeros(n^2)'] |> vecs
at = Atoms(:X, xvecs)
# set_pbc!(at, (true, true, true))
set_constraint!(at, FixedCell2D(at))
set_calculator!(at, lennardjones())
V = lennardjones()

# info("Profiling the standard LJ calculator:")
# print("Neighbourlist: "); @btime neighbourlist($at, $(cutoff(V)))
# print("energy: "); @btime energy($V, $at);
# print("forces: "); @btime forces($V, $at);

info("Profile the Experimental calculator")
ljf = Experiment.FastLJ()
energy(ljf, at) # warmup
at = Atoms(:X, xvecs) # new at to clear the data
print("First energy: "); @time energy(ljf, at) # first assembly requires nlist
print("Second energy: "); @btime energy(ljf, at) # second assembly uses stored nlist

# allocate storage for the forces
F = zeros(JVecF, length(at))

info("Profiling energy, forces, forces! without assembly:")
print("energy: "); @btime energy($ljf, $at)
print("forces: "); @btime forces($ljf, $at)
print("forces!: "); @btime Experiment.forces!($F, $ljf, $at);

# pseudo-evolution: should observe occasional nlist assembly
# for n = 1:50
#     rattle!(at, 0.1)
#     @time (energy(ljf, at), forces(ljf, at));
# end


text:
include("/home/user/Documents/HessianCalc/FFCalc.jl")

function parser(filename::String, Natoms::Int64)
    if filename[end-3:end] == "hess"
        try
            open(filename, "r") do io
                Nfree = Natoms*3
                hessian = zeros(Float64, Nfree, Nfree)
                re5 = r"[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)\s*[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)\s*[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)\s*[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)\s*[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)"
                re4 = r"[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)\s*[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)\s*[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)\s*[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)"
                re3 = r"[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)\s*[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)\s*[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)"
                re2 = r"[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)\s*[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)"
                re1 = r"[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)"
                flag = false
                i = 1
                j = 1
                for ln in eachline(io)
                    if ln === "\$hessian"
                        flag = true
                    elseif ln === "\$vibrational_frequencies"
                        break
                    elseif flag === true & occursin(re1, ln) === true
                        if length(rsplit(ln)) === 6
                            hess_str = match(re5, ln)
                            hess_val = rsplit(hess_str.match)
                            hess1 = parse(Float64, hess_val[1]); hess2 = parse(Float64, hess_val[2]);
                            hess3 = parse(Float64, hess_val[3]); hess4 = parse(Float64, hess_val[4]);
                            hess5 = parse(Float64, hess_val[5]);
                            hessian[i, j] = hess1; hessian[i, j + 1] = hess2;
                            hessian[i, j + 2] = hess3; hessian[i, j + 3] = hess4;
                            hessian[i, j + 4] = hess5;

                        elseif length(rsplit(ln)) === 5
                            hess_str = match(re4, ln)
                            hess_val = rsplit(hess_str.match)
                            hess1 = parse(Float64, hess_val[1]); hess2 = parse(Float64, hess_val[2]);
                            hess3 = parse(Float64, hess_val[3]); hess4 = parse(Float64, hess_val[4]);
                            hessian[i, j] = hess1; hessian[i, j + 1] = hess2;
                            hessian[i, j + 2] = hess3; hessian[i, j + 3] = hess4;

                        elseif length(rsplit(ln)) === 4
                            hess_str = match(re3, ln)
                            hess_val = rsplit(hess_str.match)
                            hess1 = parse(Float64, hess_val[1]); hess2 = parse(Float64, hess_val[2]);
                            hess3 = parse(Float64, hess_val[3]);
                            hessian[i, j] = hess1; hessian[i, j + 1] = hess2;
                            hessian[i, j + 2] = hess3;

                        elseif length(rsplit(ln)) === 3
                            hess_str = match(re2, ln)
                            hess_val = rsplit(hess_str.match)
                            hess1 = parse(Float64, hess_val[1]); hess2 = parse(Float64, hess_val[2]);
                            hessian[i, j] = hess1; hessian[i, j + 1] = hess2;

                        elseif length(rsplit(ln)) === 2
                            hess_str = match(re1, ln)
                            hess_val = rsplit(hess_str.match)
                            hess1 = parse(Float64, hess_val[1]);
                            hessian[i, j] = hess1;

                        else
                            continue
                        end

                        j = j + div(i, Nfree) * 5
                        i = i%(Nfree) + 1
                    else
                        continue
                    end
                end
                return(hessian)
            end
        catch
            @warn "File not founded or wrong parsing"
        end

    elseif filename[end-2:end] == "xyz"
        try
            open(filename, "r") do io
                xyz = zeros(Float64, Natoms, 3)
                re = r"([+-]?[0-9]*[.])?[0-9]+\s*[+-]?([0-9]*[.])?[0-9]+\s*[+-]?([0-9]*[.])?[0-9]+(\s*)"
                i = 1
                for ln in eachline(io)
                    if occursin(re, ln) === true
                        m = match(re, ln)
                        coord = rsplit(m.match)
                        xyz[i, 1] = parse(Float64, coord[1])
                        xyz[i, 2] = parse(Float64, coord[2])
                        xyz[i, 3] = parse(Float64, coord[3])
                        i += 1
                    else
                        continue
                    end
                end
                return(xyz)
            end
        catch
            @warn "File not founded or wrong parsing"
        end
    else
        println("Warning!!! Wrong file extension")
    end
end

xyz = parser("/home/user/Documents/HessianCalc/benzene.xyz", 12);
hessian = parser("/home/user/Documents/HessianCalc/benzene.hess", 12);
hessian = round.(hessian; digits=2)
xyz = round.(xyz; digits=2)
system = Hess_data(hessian, xyz)
bonds = [2 1; 3 2; 4 3; 5 4; 6 1; 7 1; 8 2; 9 3; 10 4; 11 5; 12 6; 6 5]
angles  = [1 2 3; 2 3 4;3 4 5;2 1 6;2 1 7;1 2 8;2 3 9;
3 4 10;4 5 11;1 6 12;6 5 11;6 1 7;5 6 12;5 4 10;1 6 5;
4 3 9;4 5 6;3 2 8;]
impropers = [12     6     1     5;
11     5     4     6;
10     4     3     5;
 9     3     2     4;
 8     2     1     3;
 7     1     2     6;]

println("**************************************")
println("*********   BONDS CONSTANTS   ********")
println("**************************************")
for i in 1:size(bonds)[1]
    println(bonds[i, 1], " ", bonds[i, 2], " bond constant:  ", bond_constant(bonds[i, 1] , bonds[i, 2] , system))
end

println("**************************************")
println("*********  ANGLES CONSTANTS   ********")
println("**************************************")
for i in 1:size(angles)[1]
    println(angles[i, 1], " ", angles[i, 2], " ",  angles[i, 3], " angle constant:  ", angle_constant(angles[i, 1] , angles[i, 2], angles[i, 3] , system))
end

println("**************************************")
println("******  IMPROPERS CONSTANTS   ********")
println("**************************************")
for i in 1:size(impropers)[1]
    println(impropers[i, 1], " ", impropers[i, 2], " ",  impropers[i, 3]," ",  impropers[i, 4],  " angle constant:  ", improper_constant(impropers[i, 1], impropers[i, 2], impropers[i, 3], impropers[i, 4] , system))
end
println("**************************************")

max_stars_repo_name: Dimitrifulgur/SeminarioFUERZA

max_issues_repo_name: Dimitrifulgur/SeminarioFUERZA

max_forks_repo_name: Dimitrifulgur/SeminarioFUERZA


text:
function initial_velocity(N::Int64;T=298.15,R=8.3145,m=1.0)
  vel = [ [0.,0.] for i in 1:N]

  for i in 1:N
    angle = 2*pi*rand()
    px = rand()
    py = rand()
    vel[i][1] = sqrt(-(2*R*T/m)*log(1-px)) * angle       # velocity in the x axis
    vel[i][2] = sqrt(-(2*R*T/m)*log(1-py)) * angle       # velocity in the y axis
  end

  return vel
end

function norm(x::Vector{Float64})
  s = 0.
  for i in 1:length(x)
    s = s + x[i]^2
  end
  sqrt(s)
end

function normVEL(x)
  s = zeros(length(x))
  for i in 1:length(x)
    s[i] = norm(x[i])
  end
  return s
end

export initial_velocity, norm, normVEL

# p(v)dv
p(v;m=1.,R=8.3145,T=298.15) = m/(R*T) * v * exp(-m*v^2/(2R*T))

# probability v < x
pv(x;m = 1., T = 298.15, R = 8.3145) = -exp(-m*x^2/(2*R*T)) + 1


text:
function prepare_spatial_data(input_data)
    # OHL
    spatial_data_ohl, costs_ohl, equipment_data_ohl = build_spatial_matrices_and_cost_data(input_data; equipment = "ohl")
    spatial_data_ugc, costs_ugc, equipment_data_ugc = build_spatial_matrices_and_cost_data(input_data; equipment = "ugc")

    spatial_data = merge_maps(spatial_data_ohl, costs_ohl, spatial_data_ugc, costs_ugc, input_data)

    spatial_data_matrices = Dict{String, Any}()
    cost_data = Dict{String, Any}()
    equipment_data = Dict{String, Any}()
    spatial_data_matrices["ohl"] = spatial_data_ohl
    spatial_data_matrices["ugc"] = spatial_data_ugc
    cost_data["ohl"] = costs_ohl
    cost_data["ugc"] = costs_ugc
    equipment_data["ohl"] = equipment_data_ohl
    equipment_data["ugc"] = equipment_data_ugc

    return spatial_data, spatial_data_matrices, cost_data, equipment_data
end

function build_spatial_matrices_and_cost_data(input_data; equipment = "ohl")
    # prepare boundaries
    xmax = size(input_data["rgb_values"]["sea"], 2) + 1
    ymax = size(input_data["rgb_values"]["sea"], 1) + 1

    x_border_min = Int(input_data["boundaries"]["xmin"])
    x_border_max = Int(input_data["boundaries"]["xmax"])
    y_border_min = Int(input_data["boundaries"]["ymin"])
    y_border_max = Int(input_data["boundaries"]["ymax"])

    x_node_min = Int(min(input_data["start_node"]["x"], input_data["end_node"]["x"]))
    x_node_max = Int(max(input_data["start_node"]["x"], input_data["end_node"]["x"]))
    y_node_min = Int(min(input_data["start_node"]["y"], input_data["end_node"]["y"]))
    y_node_max = Int(max(input_data["start_node"]["y"], input_data["end_node"]["y"]))

    delta_x = Int(max(input_data["resolution_factor"],(x_node_max-x_node_min)/xmax))
    delta_y = Int(max(input_data["resolution_factor"],(y_node_max-y_node_min)/ymax))

    # make grid

    x = collect(x_node_min : delta_x : x_node_max)
    y = collect(y_node_min : delta_y : y_node_max)

    if x[end] != x_node_max
        x = [x; x_node_max]
    end
    if y[end] != y_node_max
        y = [y; y_node_max]
    end

    x_beg = collect(x_border_min : delta_x : x[1])
    if x_beg[end] != x[1]
        x_beg = [x; x[1]]
    end
    if length(x) == 1
        x_beg[end] = []
    end
    y_beg = collect(y_border_min : delta_y : y[1])
    if y_beg[end] != y[1]
        y_beg = [y; y[1]]
    end
    if length(y) == 1
        y_beg[end]= []
    end

    x = [x_beg; x[2:end-1]; collect(x[end] : delta_x : x_border_max)]
    y = [y_beg; y[2:end-1]; collect(y[end] : delta_y : y_border_max)]

    X = x' .* ones(size(y, 1))
    Y = (ones(size(x, 1)))' .*  y

    node_weights = Dict{String, Any}()
    for (area_idx, area) in input_data["rgb_values"]
        node_weights[area_idx] = assign_rgb_values(X, Y, x, y, input_data, area_idx)
    end

    onshore_nodes = ones(size(node_weights["sea"],1),1)
    onshore_nodes = 1 .- node_weights["sea"][:,4]

    nodes = node_weights["sea"][:, 1:3]

    equipment_ = join([input_data["technology"],"_",equipment])
    node_weights = assign_node_weights(node_weights, input_data, equipment_)

    segments, segment_weights, segment_km, onshore_flag = create_segments_and_segment_weights(nodes, node_weights, onshore_nodes, input_data, X, Y, x, y, delta_x, delta_y)

    costs, equipment_details = calculate_costs_and_losses(input_data; equipment = equipment)

    spatial_data = Dict{String, Any}()
    spatial_data["nodes"] = nodes
    spatial_data["node_weights"] = node_weights
    spatial_data["onshore_nodes"] = onshore_nodes
    spatial_data["segments"] = segments
    spatial_data["segment_weights"] = segment_weights
    spatial_data["segment_km"] = segment_km
    spatial_data["onshore_flag"] = onshore_flag
    return spatial_data, costs, equipment_details
end

function assign_rgb_values(X, Y, x, y, input_data, area)
    nodes_w = zeros(size(x, 1) * size(y, 1), 4)
    for idx = 1 : size(X, 2)
        nodes_w[((idx-1) * size(X,1)+1) : idx * size(X,1) ,1] = ((idx-1) * size(X,1)+1) :idx * size(X,1)
        nodes_w[((idx-1) * size(X,1)+1) : idx * size(X,1), 2] = X[:, idx]
        nodes_w[((idx-1) * size(X,1)+1) : idx * size(X,1), 3] = Y[:, idx]
        nodes_w[((idx-1) * size(X,1)+1) : idx * size(X,1), 4] = input_data["rgb_values"][area][y, Int.(x[idx])]
    end
    return nodes_w
end

function assign_node_weights(node_weights, input_data, equipment)

    for (area, value) in input_data["spatial_weights"][equipment]
        node_weights[area][:, 4] = node_weights[area][:, 4] * input_data["spatial_weights"][equipment][area]
    end

    node_weights_ = zeros(size(node_weights["sea"]))
    node_weights_[:, 1:3] = node_weights["sea"][:, 1:3]
    for idx = 1:size(node_weights_, 1)
        for (area, value) in node_weights
            if node_weights_[idx, 4] == 0 && node_weights[area][idx, 4] !=0
                node_weights_[idx, 4] = node_weights[area][idx, 4]
            end
            if node_weights_[idx, 4] != 0 && node_weights[area][idx, 4] !=0
                if input_data["overlapping_area_weight"] == "minimum"
                    node_weights_[idx, 4] = min(node_weights[area][idx, 4], node_weights_[idx, 4])
                else
                    node_weights_[idx, 4] = (node_weights[area][idx, 4] + node_weights_[idx, 4]) / 2
                end
            end
        end
        node_weights_[idx, 4] = max(1, node_weights_[idx, 4])
    end
    return node_weights_
end

function create_segments_and_segment_weights(nodes, node_weights, onshore_nodes, input_data, X, Y, x, y, delta_x, delta_y)

    segments = zeros((length(x)-1)*length(y)+(length(y)-1)*length(x)+2*((length(x)-1)*(length(y)-1)),3);
    segments[:,1]=(1:size(segments,1))'
    segment_weights=ones(size(segments,1),1)
    onshore_flag=ones(size(segments,1),1)
    segment_km=ones(size(segment_weights))

    # Construct all segments: Up, to the side, and diagonal up.
    for idx = 1 : size(X, 2) - 1
        # PART1
        segments[(4*(idx-1)*size(X,1)+1-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2] = nodes[((idx-1)*size(X,1)+1):idx*size(X,1),1]
        segments[(4*(idx-1)*size(X,1)+1-3*(idx-1)):4:4*idx*size(X,1)-3*idx,3] = segments[(4*(idx-1)*size(X,1)+1-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2] .+ 1
        if input_data["overlapping_area_weight"] == "minimum"
            segment_weights[(4*(idx-1)*size(X,1)+1-3*(idx-1)):4:4*idx*size(X,1)-3*idx] = minimum(node_weights[((idx-1)*size(X,1)+1):idx*size(X,1), 4],node_weights[segments[(4*(idx-1)*size(X,1)+1-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2+1],4])
        elseif input_data["overlapping_area_weight"] == "average"
            w1 = node_weights[((idx-1)*size(X,1)+1):idx*size(X,1),4]
            w2 = node_weights[Int.(segments[(4*(idx-1)*size(X,1)+1-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2] .+ 1), 4]
            segment_weights[(4*(idx-1)*size(X,1)+1-3*(idx-1)):4:4*idx*size(X,1)-3*idx] = (w1 + w2) ./ 2
        else
            print("Overlapping area weight must be average or minimum")
        end
        idx_ = Int.(segments[(4*(idx-1)*size(X,1)+1-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2].+1)
        idx__ = ((idx-1)*size(X,1)+1):idx*size(X,1)
        onshore_flag[(4*(idx-1)*size(X,1)+1-3*(idx-1)):4:4*idx*size(X,1)-3*idx] = minimum([onshore_nodes[idx__],onshore_nodes[idx_]])
        segments[4*idx*size(X,1)-3*idx,3] = segments[4*idx*size(X,1)-3*idx,2]+size(X,1)
        idx1 = (4*(idx-1)*size(X,1)+1-3*(idx-1)):4:4*idx*size(X,1)-3*idx
        segment_km[idx1] = sqrt.((nodes[Int.(segments[idx1,2]),2]-nodes[Int.(segments[idx1,3]),2]).^2+(nodes[Int.(segments[idx1,2]),3]-nodes[Int.(segments[idx1,3]),3]).^2)
        # PART 2
        segments[(4*(idx-1)*size(X,1)+2-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2] = nodes[((idx-1)*size(X,1)+1):idx*size(X,1)-1,1]
        segments[(4*(idx-1)*size(X,1)+2-3*(idx-1)):4:4*idx*size(X,1)-3*idx,3] = segments[(4*(idx-1)*size(X,1)+2-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2] .+ size(X,1)
        if input_data["overlapping_area_weight"] == "minimum"
            segment_weights[(4*(idx-1)*size(X,1)+2-3*(idx-1)):4:4*idx*size(X,1)-3*idx] = mininimum(node_weights[((idx-1)*size(X,1)+1):idx*size(X,1)-1, 4],node_weights[segments[(4*(idx-1)*size(X,1)+2-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2]+size(X,1), 4])
        elseif input_data["overlapping_area_weight"] == "average"
            w1 = node_weights[((idx-1)*size(X,1)+1):idx*size(X,1)-1, 4]
            w2 = node_weights[Int.(segments[(4*(idx-1)*size(X,1)+2-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2] .+ size(X,1)), 4]
            segment_weights[(4*(idx-1)*size(X,1)+2-3*(idx-1)):4:4*idx*size(X,1)-3*idx] = (w1 + w2) ./ 2
        else
            print("Overlapping area weight must be average or minimum")
        end
        idx__ = Int.(((idx-1)*size(X,1)+1):idx*size(X,1) .- 1)
        idx_ = Int.(segments[(4*(idx-1)*size(X,1)+2-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2] .+ size(X,1))
        onshore_flag[(4*(idx-1)*size(X,1)+2-3*(idx-1)):4:4*idx*size(X,1)-3*idx] = minimum([onshore_nodes[idx__],onshore_nodes[idx_]])
        idx1 = Int.((4*(idx-1)*size(X,1)+2-3*(idx-1)):4:4*idx*size(X,1)-3*idx)
        segment_km[idx1] = sqrt.((nodes[Int.(segments[idx1,2]),2]-nodes[Int.(segments[idx1,3]),2]).^2+(nodes[Int.(segments[idx1,2]),3]-nodes[Int.(segments[idx1,3]),3]).^2)
        #PART 3
        segments[(4*(idx-1)*size(X,1)+3-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2] = nodes[((idx-1)*size(X,1)+1):idx*size(X,1)-1,1]
        segments[(4*(idx-1)*size(X,1)+3-3*(idx-1)):4:4*idx*size(X,1)-3*idx,3] = segments[(4*(idx-1)*size(X,1)+3-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2] .+ size(X,1) .+ 1
        if input_data["overlapping_area_weight"] == "minimum"
            segment_weights[(4*(idx-1)*size(X,1)+3-3*(idx-1)):4:4*idx*size(X,1)-3*idx] = minimum(node_weights[((idx-1)*size(X,1)+1):idx*size(X,1)-1, 4], node_weights[segments[(4*(idx-1)*size(X,1)+3-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2] + size(X,1)+1, 4])
        elseif input_data["overlapping_area_weight"] == "average"
            w1 = node_weights[((idx-1)*size(X,1)+1):idx*size(X,1)-1, 4]
            w2 = node_weights[Int.(segments[(4*(idx-1)*size(X,1)+3-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2] .+ size(X,1) .+ 1), 4]
            segment_weights[(4*(idx-1)*size(X,1)+3-3*(idx-1)):4:4*idx*size(X,1)-3*idx] = (w1 + w2) ./ 2
        else
            print("Overlapping area weight must be average or minimum")
        end
        idx__ = Int.(((idx-1)*size(X,1)+1):idx*size(X,1)-1)
        idx_ = Int.(segments[(4*(idx-1)*size(X,1)+3-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2] .+ size(X,1) .+ 1)
        onshore_flag[(4*(idx-1)*size(X,1)+3-3*(idx-1)):4:4*idx*size(X,1)-3*idx] = minimum([onshore_nodes[idx__],onshore_nodes[idx_]])
        idx1 = Int.((4*(idx-1)*size(X,1)+3-3*(idx-1)):4:4*idx*size(X,1)-3*idx)
        segment_km[idx1]=sqrt.((nodes[Int.(segments[idx1,2]),2]-nodes[Int.(segments[idx1,3]),2]).^2+(nodes[Int.(segments[idx1,2]),3]-nodes[Int.(segments[idx1,3]),3]).^2)
        #PART 4
        segments[(4*(idx-1)*size(X,1)+4-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2] = nodes[((idx-1)*size(X,1)+1)+1:idx*size(X,1),1]
        segments[(4*(idx-1)*size(X,1)+4-3*(idx-1)):4:4*idx*size(X,1)-3*idx,3] = segments[(4*(idx-1)*size(X,1)+2-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2] .+ size(X,1)
        if input_data["overlapping_area_weight"] == "minimum"
            segment_weights[(4*(idx-1)*size(X,1)+4-3*(idx-1)):4:4*idx*size(X,1)-3*idx] = minimum(node_weights[((idx-1)*size(X,1)+1)+1:idx*size(X,1), 4],node_weights[segments[(4*(idx-1)*size(X,1)+2-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2]+size(X,1), 4])
        elseif input_data["overlapping_area_weight"] == "average"
            w1 = node_weights[((idx-1)*size(X,1)+1)+1:idx*size(X,1), 4]
            w2 = node_weights[Int.(segments[(4*(idx-1)*size(X,1)+2-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2] .+ size(X,1)), 4]
            segment_weights[(4*(idx-1)*size(X,1)+4-3*(idx-1)):4:4*idx*size(X,1)-3*idx] = (w1 + w2) ./ 2
        else
            print("Overlapping area weight must be average or minimum")
        end
        idx__ = ((idx-1)*size(X,1)+1)+1:idx*size(X,1)
        idx_ = Int.(segments[(4*(idx-1)*size(X,1)+2-3*(idx-1)):4:4*idx*size(X,1)-3*idx,2] .+ size(X,1))
        onshore_flag[(4*(idx-1)*size(X,1)+4-3*(idx-1)):4:4*idx*size(X,1)-3*idx] = minimum([onshore_nodes[idx__],onshore_nodes[idx_]])
        idx1 = Int.((4*(idx-1)*size(X,1)+4-3*(idx-1)):4:4*idx*size(X,1)-3*idx)
        segment_km[idx1]=sqrt.((nodes[Int.(segments[idx1,2]),2]-nodes[Int.(segments[idx1,3]),2]).^2+(nodes[Int.(segments[idx1,2]),3]-nodes[Int.(segments[idx1,3]),3]).^2)
    end

    segments[end-size(X,1)+2:end,2] = nodes[end-size(X,1)+1:end-1,1]
    segments[end-size(X,1)+2:end,3] = segments[end-size(X,1)+2:end,2] .+ 1
    if input_data["overlapping_area_weight"] == "minimum"
        segment_weights[end-size(X,1)+2:end] = minimum(node_weights[end-size(X,1)+1:end-1, 4], node_weights[Int.(segments[end-size(X,1)+2:end, 2] .+ 1) , 4])
    elseif input_data["overlapping_area_weight"] == "average"
        w1 = node_weights[end-size(X,1)+1:end-1, 4]
        w2 = node_weights[Int.(segments[end-size(X,1)+2:end,2] .+ 1), 4]
        segment_weights[end-size(X,1)+2:end] = (w1 + w2) ./ 2
    else
        print("Overlapping area weight must be average or minimum")
    end
    onshore_flag[end-size(X,1)+2:end] = onshore_nodes[Int.(segments[end-size(X,1)+2:end,2] .+ 1)] .* onshore_nodes[end - size(X,1) + 1 : end-1]
    segment_km[end - size(X,1) + 2 : end] .= delta_y * input_data["resolution_factor"]

    return segments, segment_weights, segment_km, onshore_flag
end

function merge_maps(spatial_data_ohl, costs_ohl, spatial_data_ugc, costs_ugc, input_data)
    costs_ohl["ohl_km"] = spatial_data_ohl["segment_km"] .* (spatial_data_ohl["segment_weights"] .* costs_ohl["installation"] .+ costs_ohl["investment"]) ./ input_data["resolution_factor"]
    costs_ugc["ugc_km"] = spatial_data_ugc["segment_km"] .* ((spatial_data_ugc["segment_weights"] .* ((costs_ugc["onshore"]["installation"] .* spatial_data_ugc["onshore_flag"]) .+ (costs_ugc["offshore"]["installation"] .* (1 .- spatial_data_ugc["onshore_flag"])))) .+ ((costs_ugc["onshore"]["investment"] .* spatial_data_ugc["onshore_flag"]) .+ (costs_ugc["offshore"]["investment"] .* (1 .- spatial_data_ugc["onshore_flag"])))) ./ input_data["resolution_factor"]

    # Update segment weights based on costs
    spatial_data_ohl["segment_weights"] = (spatial_data_ohl["segment_weights"] .* costs_ohl["installation"] .+ costs_ohl["investment"]) .* input_data["resolution_factor"]
    spatial_data_ugc["segment_weights"] = ((spatial_data_ugc["segment_weights"] .* ((costs_ugc["onshore"]["installation"] .* spatial_data_ugc["onshore_flag"]) .+ (costs_ugc["offshore"]["installation"] .* (1 .- spatial_data_ugc["onshore_flag"])))) .+ ((costs_ugc["onshore"]["investment"] .* spatial_data_ugc["onshore_flag"]) .+ (costs_ugc["offshore"]["investment"] .* (1 .- spatial_data_ugc["onshore_flag"])))) .* input_data["resolution_factor"]

    # # Connect OHL and UGC maps
    # # Step 1) Shift node indexes for the UGC map with the number of nodes for OHL
    spatial_data_ugc["segments"][:, 1] = spatial_data_ugc["segments"][:, 1] .+ size(spatial_data_ohl["segments"], 1)
    spatial_data_ugc["segments"][:, 2:3] = spatial_data_ugc["segments"][:, 2:3] .+ size(spatial_data_ohl["nodes"], 1)
    spatial_data_ugc["nodes"][:, 1] = spatial_data_ohl["nodes"][:,1] .+ size(spatial_data_ohl["nodes"], 1)

    spatial_data = Dict{String, Any}()
    spatial_data["segments"] = [spatial_data_ohl["segments"]; spatial_data_ugc["segments"]]
    spatial_data["nodes"] = [spatial_data_ohl["nodes"]; spatial_data_ugc["nodes"]]
    spatial_data["onshore_nodes"] = [spatial_data_ohl["onshore_nodes"]; spatial_data_ugc["onshore_nodes"]]
    spatial_data["segment_weights"] = [spatial_data_ohl["segment_weights"]; spatial_data_ugc["segment_weights"]]
    spatial_data["onshore_flag"] = [spatial_data_ohl["onshore_flag"]; spatial_data_ugc["onshore_flag"]]

    # # Step 2) Connect OHL and UGC maps at each X, Y point (node) via a new segment
    number_of_additional_segments = size(spatial_data_ohl["nodes"], 1)

    segment_idx = 1 : (spatial_data["segments"][end, 1] + number_of_additional_segments)
    segment_from_node =  [spatial_data["segments"][:, 2]; spatial_data_ohl["nodes"][:, 1]]
    segment_to_node = [spatial_data["segments"][:, 3]; spatial_data_ugc["nodes"][:, 1]]

    spatial_data["segments"] = [segment_idx segment_from_node segment_to_node]
    spatial_data["segment_weights"] = [spatial_data["segment_weights"]; zeros(number_of_additional_segments, 1)]
    spatial_data["onshore_flag"] = [spatial_data["onshore_flag"]; ones(number_of_additional_segments, 1)]
    spatial_data["ohl_ugc_transition"] = [zeros(2 * size(spatial_data["segments"], 1)); ones(number_of_additional_segments, 1)]

    return spatial_data
end


text:
"""
    symbol(plugin)

Return the per-PluginStack unique Symbol of this plugin if it exports a "late-bind" runtime API to other plugins.
"""
symbol(plugin::Plugin) = :nothing

"""
    setup!(plugin, deps, args...)

Initialize the plugin with the given dependencies and arguments (e.g. shared state).

This lifecycle hook will be called when the application loads a plugin. Plugins.jl does not (yet) helps with this,
application developers should do it manually, right after the PluginStack was created, before the hook_cache() call.
"""
setup!(plugin::Plugin, args...) = nothing

"""
    shutdown!(plugin, args...)

Shut down the plugin.

This lifecycle hook will be called when the application unloads a plugin,
e.g. before the application exits.
Plugins.jl does not (yet) helps with this, application developers should do it manually.
"""
shutdown!(plugin::Plugin, args...) = nothing

"""
    customfield(plugin::Plugin, abstract_type::Type, args...) = nothing

Provide field specifications to plugin-assembled types.

Using this lifecycle hook the system can define custom plugin-assembled types
(typically structs) based on field specifications provided by plugins.
E.g. an error type can be extended with debug information.

A plugin can provide zero or one field to every assembled type.
To provide a field, return a [`FieldSpec`](@ref).

The assembled type will be a subtype of `abstract_type`. To allow differently
configured systems to run in the same Julia session, new types may be assembled
for every instance of the system.

!!! warning "Metaprogramming may make you unhappy"
    Although plugin-assmebled types are designed to help doing metaprogramming in a controlled
    fashion, it is usually better to use non-meta solutions instead. E.g. Store plugin state
    inside the plugin, collect data from multiple plugins using lifecycle hooks, etc.
"""
customfield(plugin::Plugin, abstract_type::Type, args...) = nothing

"""
    request_stage(plugin::Plugin, args...)::Stage

Request a new stage iteration by returning a Stage representing it.

This lifecycle hook will be called repeatedly to ask plugins
their wish to stage. If a plugin returns a `Stage` instance
and the request is accepted, the stage will start immediately.

If more than one plugins ask for staging, their request will
be merged if possible and only one stage will run. If the stages
are incompatible, meaning that different sets of plugins handle the
`prepeare` hook of the stages, then only a compatible subset
of them will run.

Plugins should continue requesting staging until their wish gets
fulfilled.
"""
request_stage(plugin::Plugin, args...) = nothing

"""
    prepare_stage(plugin::Plugin, stage::Stage)

Lifecycle hook to prepare the plugin for the next stage.
If the stage is an `EvalStage` and the plugin needs to evaluate
code, this is the point to do it.
"""
prepare_stage(plugin::Plugin, stage, args...) = nothing

"""
    enter_stage(plugin::Plugin, stage::Stage, args...)

Lifecycle hook marking the start of the next stage.

Types are reasssembled at this point. If `stage isa EvalStage`,
then the world is already updated. (execution reached toplevel)

"""
enter_stage(plugin::Plugin, stage, args...) = nothing

"""
    leave_stage(plugin::Plugin, stage::Stage, nextstage::Stage, args...)

Lifecycle hook marking the end of a stage.
"""
leave_stage(plugin::Plugin, stage, args...) = nothing

setup_hook! = create_lifecyclehook(setup!)
shutdown_hook! = create_lifecyclehook(shutdown!)
customfield_hook = create_lifecyclehook(customfield)
request_stage_hook = create_lifecyclehook(request_stage)
prepare_stage_hook = create_lifecyclehook(prepare_stage)
enter_stage_hook = create_lifecyclehook(enter_stage)
leave_stage_hook = create_lifecyclehook(leave_stage)

setup!(stack::PluginStack, sharedstate) = setup_hook!(stack, sharedstate)
shutdown!(stack::PluginStack, sharedstate) = shutdown_hook!(stack, sharedstate)
customfields(stack::PluginStack, abstract_type::Type) = customfield_hook(stack, abstract_type)


text:
"""
    load(path::AbstractString)

Load STL file to list of triangles.
"""
function load(path::AbstractString)
    open(path) do stl
        skip(stl, 80)  # skip header
        trianglecount = read(stl, UInt32)
        ref = Ref{Triangle}()
        triangles = map(1:trianglecount) do i
            read!(stl, ref)
            skip(stl, 2)  # skip attribute byte count
            ref[]
        end
        @assert eof(stl)
        return triangles
    end
end


text:
# ============================================================
# Thermodynamics
# ============================================================

"""
$(SIGNATURES)

Calculate heat capacity ratio (monatomic gas)
"""
function heat_capacity_ratio(K, D::T) where {T<:Integer}
    γ = begin
        if D == 1
            (K + 3.0) / (K + 1.0)
        elseif D == 2
            (K + 4.0) / (K + 2.0)
        elseif D == 3
            (K + 5.0) / (K + 3.0)
        end
    end

    return γ
end

"""
$(SIGNATURES)

Calculate heat capacity ratio (diatomic gas)
"""
function heat_capacity_ratio(K, Nr, D::T) where {T<:Integer}
    γ = begin
        if D == 1
            (K + 3.0 + Nr) / (K + 1.0 + Nr)
        elseif D == 2
            (K + 4.0 + Nr) / (K + 2.0 + Nr)
        elseif D == 3
            (K + 5.0 + Nr) / (K + 3.0 + Nr)
        end
    end

    return γ
end

"""
$(SIGNATURES)

Calculate speed of sound
"""
sound_speed(λ::Real, γ::Real) = (0.5 * γ / λ)^0.5

"""
$(SIGNATURES)
"""
sound_speed(prim::AV, γ) = sound_speed(prim[end], γ)

"""
$(TYPEDSIGNATURES)

Calculate sound speed in mixture
"""
function sound_speed(prim::AM, γ)
    c = similar(prim, axes(prim, 2))
    for j in eachindex(c)
        c[j] = sound_speed(prim[end, j], γ)
    end

    return maximum(c)
end


text:
export default_type, default_library, similar_library, library, dimension_names

"""
    default_type(d::FullDim, ::Type{T}) where {T}

Returns the default polyhedron type for `d`-dimensional polyhedron of coefficient type `T`.
"""
function default_type end

function default_type(d::StaticArrays.Size{N}, ::Type{T}) where {N, T}
    return DefaultPolyhedron{T, Intersection{T, StaticArrays.SVector{N[1], T}, typeof(d)}, Hull{T, StaticArrays.SVector{N[1], T}, typeof(d)}}
end
function default_type(d::StaticArrays.Size{(1,)}, ::Type{T}) where T
    return Interval{T, StaticArrays.SVector{1, T}, typeof(d)}
end
function default_type(d::Integer, ::Type{T}) where T
    if isone(d)
        return Interval{T, Vector{T}, typeof(d)}
    else
        return DefaultPolyhedron{T, Intersection{T, Vector{T}, typeof(d)}, Hull{T, Vector{T}, typeof(d)}}
    end
end

"""
    default_library(d::FullDim, ::Type{T}) where {T}

Returns the default polyhedral library for `d`-dimensional polyhedron of coefficient type `T`.

## Examples

To obtain the default library for 2-dimensional polyhedra of eltype `Float64`,
do `default_library(2, Float64)`.

Given an `StaticArrays.SVector` `v`, to obtain a default library for points
of the type of `v` in a type stable way, do
`default_library(Polyhedra.FullDim(v), eltype(v))`.
"""
function default_library end

_default_type(::Type{T}) where T = T
# See https://github.com/JuliaPolyhedra/Polyhedra.jl/issues/35
_default_type(::Type{AbstractFloat}) = Float64
default_library(::StaticArrays.Size, T::Type) = DefaultLibrary{_default_type(T)}()
default_library(::StaticArrays.Size{(1,)}, T::Type) = IntervalLibrary{_default_type(T)}()
function default_library(d::Integer, ::Type{T}) where T
    if isone(d)
        return IntervalLibrary{_default_type(T)}()
    else
        return DefaultLibrary{_default_type(T)}()
    end
end

"""
    similar_library(lib::Library, d::FullDim, T::Type)

Returns a library that supports polyhedra of full dimension `T` with coefficient type `T`. If `lib` does not support it, this commonly calls `default_library(d, T)`.
"""
function similar_library end

# Shortcuts
similar_library(p::Polyhedron, d::FullDim, ::Type{T}) where T = similar_library(library(p), d, T)
similar_library(p::Polyhedron, ::Type{T}) where T = similar_library(library(p), T)
similar_library(p::Polyhedron{T}, d::FullDim) where T = similar_library(p, d, T)

"""
    library(p::Polyhedron)

Returns the library used by `p`.
"""
function library end

"""
    default_solver(p::Rep)

Returns a default linear programming solver for the polyhedron `p` (e.g. CDD has an internal solver which is used by default).
"""
default_solver(p::Rep; kws...) = nothing
function default_solver(p::Rep, ps::Rep...; kws...)
    s = default_solver(p; kws...)
    if s === nothing
        return default_solver(ps...; kws...)
    else
        return s
    end
end

"""
    linear_objective_solver(p::Rep, solver=default_solver(p))

Return the solver to use for optimizing a linear objective over the polyhedron `p`, i.e.
```julia
model = Model(solver)
x = @variable(model, [1:fulldim(p)])
@constraint(model, x in p)
@objective(model, c ⋅ x)
```
for some vector `c`.

By default, if the V-representation of `p` has been computed, it returns
`VRepOptimizer()`, otherwise, it returns `solver`.

If the problem has constraints different to `x in p`, use `default_solver(p)` instead
as the fact that the V-representation of `p` has been computed does not help.
"""
function linear_objective_solver(p::Rep{T}, solver=default_solver(p)) where T
    if vrepiscomputed(p)
        return VRepOptimizer{T}
    else
        return solver
    end
end

linear_objective_solver(v::VRepresentation{T}, solver=default_solver(v)) where {T} = VRepOptimizer{T}
linear_objective_solver(h::HRepresentation, solver=default_solver(h)) = solver

_promote_reptype(P::Type{<:HRep}, ::Type{<:HRep}) = P
_promote_reptype(P::Type{<:VRep}, ::Type{<:VRep}) = P
# Breaks ambiguity with above two methods
_promote_reptype(P::Type{<:Polyhedron}, ::Type{<:Polyhedron}) = P

function promote_reptype(P1::Type{<:Rep}, P2::Type{<:Rep})
    return _promote_reptype(P1, P2)
end
promote_reptype(P1::Type{<:Rep}, P2::Type{<:Rep}, P::Type{<:Rep}...) = promote_reptype(promote_reptype(P1, P2), P...)
promote_reptype(P::Type{<:Rep}) = P

# whether Rep has the constructor Rep(::It...; solver=...)
supportssolver(::Type{<:Rep}) = false
supports_names(::Type{<:Rep}) = false
dimension_names(::Rep) = nothing

function constructpolyhedron(RepT::Type{<:Rep{T}}, d::FullDim, p::Tuple{Vararg{Rep}}, it::It{T}...; dimension_map::Function = i -> (1, i)) where T
    solver = names = nothing
    if supportssolver(RepT)
        solver = default_solver(p...)
        if solver !== nothing
        end
    end
    if supports_names(RepT)
        names = nothing
        sub_names = dimension_names.(p)
        for i in 1:fulldim(d)
            j = dimension_map(i)
            if j !== nothing && sub_names[j[1]] !== nothing
                if names === nothing
                    names = fill("", fulldim(d))
                end
                names[i] = sub_names[j[1]][j[2]]
            end
        end
    end
    if names !== nothing
        if solver !== nothing
            return RepT(d, it...; solver=solver, dimension_names=names)::RepT
        else
            return RepT(d, it...; dimension_names=names)::RepT
        end
    else
        if solver !== nothing
            return RepT(d, it...; solver=solver)::RepT
        else
            # FIXME without this type annotation even convexhull(::PointsHull{2,Int64,Array{Int64,1}}, ::PointsHull{2,Int64,Array{Int64,1}}) is not type stable, why ?
            return RepT(d, it...)::RepT
        end
    end
end

function default_similar(p::Tuple{Vararg{Rep}}, d::FullDim, ::Type{T}, it::It{T}...; kws...) where T
    # Some types in p may not support `d` or `T` so we call `similar_type` after `promote_reptype`
    RepT = similar_type(promote_reptype(typeof.(p)...), d, T)
    return constructpolyhedron(RepT, d, p, it...; kws...)
end

"""
    similar(p::Tuple{Vararg{Polyhedra.Rep}}, d::Polyhedra.FullDim, ::Type{T}, it::Polyhedra.It{T}...)

Creates a representation with a type similar to `p` of a polyhedron of full dimension `d`, element type `T` and initialize it with the iterators `it`.
The type of the result will be chosen closer to the type of `p[1]`.
"""
Base.similar(p::Tuple{Vararg{Rep}}, d::FullDim, ::Type{T}, it::It{T}...; kws...) where {T} = default_similar(p, d, T, it...; kws...)
function promote_coefficient_type(p::Tuple{Vararg{Rep}})
    return promote_type(coefficient_type.(p)...)
end
function Base.similar(p::Tuple{Vararg{Rep}}, d::FullDim, it::It...; kws...)
    T = promote_coefficient_type(p)
    return similar(p, d, T, it...; kws...)
end
function Base.similar(p::Tuple{Vararg{Rep}}, it::It...; kws...)
    return similar(p, FullDim(p[1]), it...; kws...)
end
Base.similar(p::Rep, args...; kws...) = similar((p,), args...; kws...)

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
using DEBmicroTrait
using Test

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

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
using ReferenceMicrocalFiles
using Test

# write your own tests here
@test 1 == 1
for (k,v) in ReferenceMicrocalFiles.dict
  @test isfile(v.filename)
  @test isfile(v.noise_filename) || v.noise_filename == ""
  @test !isempty(v.note)
  @test !isempty(v.line_energies_ev)
end


text:
function MultidimensionalCounter(arg0::Vector{jint})
    return MultidimensionalCounter((Vector{jint},), arg0)
end

function get_count(obj::MultidimensionalCounter, arg0::Vector{jint})
    return jcall(obj, "getCount", jint, (Vector{jint},), arg0)
end

function get_counts(obj::MultidimensionalCounter, arg0::jint)
    return jcall(obj, "getCounts", Vector{jint}, (jint,), arg0)
end

function get_dimension(obj::MultidimensionalCounter)
    return jcall(obj, "getDimension", jint, ())
end

function get_size(obj::MultidimensionalCounter)
    return jcall(obj, "getSize", jint, ())
end

function get_sizes(obj::MultidimensionalCounter)
    return jcall(obj, "getSizes", Vector{jint}, ())
end

function iterator(obj::MultidimensionalCounter)
    return jcall(obj, "iterator", MultidimensionalCounter_JIterator, ())
end

function to_string(obj::MultidimensionalCounter)
    return jcall(obj, "toString", JString, ())
end